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ABSTRACT 


A new method has been devised to determine the 
spherical harmonic coefficients of the lunar gravity field. 

This method consists of a two-step data reduction and estima- 
tion process. In the first step, a weighted least-squares 
empirical orbit determination scheme is applied to Doppler 
tracking data from lunar orbits to estimate long-period Kepler 
elements and rates. Each of the Kepler elements is represented 
by an independent function of time. The long-period perturbing 
effects of the earth, sun, and solar radiation are explicitly 
modeled in this scheme. Kepler element variations estimated by 
this empirical processor are ascribed to the non-central lunar 
gravitation features. Doppler data are reduced in this manner 
for as many orbits as are available. In the second step, the 
Kepler element rates are used as input to a second least-squares 
processor that estimates lunar gravity coefficients using the 
long-period Lagrange perturbation equations . 

Pseudo Doppler data have been generated simulating two 


different lunar orbits. This analysis included the perturbing 
effects of triaxial lunar harmonics, the earth, sun, and solar 
radiation pressure. Orbit determinations were performed on these 
data and long-period orbital elements obtained. The Kepler 
element rates from these solutions were used to recover triaxial 
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lunar gravity coefficients. Overall results of this con- 
trolled experiment show that lunar gravity coefficients can be 
accurately determined and that the method is dynamically 
compatible with long-period perturbation theory. 

This selenodesy method has been applied to Doppler 
data from the Lunar Orbiter I, II, III, and V missions. One 
hundred ninety-nine sets of Kepler element rates are obtained 
for gravity field determination. A lunar field of degree and 
order four is derived from these rates. Equipotential surfaces 
from this gravity field show the lunar mass distribution to be 
that of a triaxial ellipsoid with three large areas of mass con- 
centration. The greatest and by far the most dominant of these 
areas is centered very near the Mare Serenitatis region and 
covers a large portion of the front side of the moon. The other 

two regions of mass concentration are located on the far side of 
the moon but do not correspond to a specific maria region. 

This gravity field has been investigated using data 
from several Apollo missions. Orbit determinations from these 
data show this field results in improved orbit predictions when 
compared to those using two other gravity fields. All solutions 
indicate the lunar gravity field models are still incomplete. 
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CHAPTER I 


INTRODUCTION 


GENERAL 

Observational measurements obtained from earth 
based tracking of satellites in close lunar orbits provide 
a unique data source for the identification of lunar gravity 
parameters. An accurate knowledge of the lunar gravity 
field is essential for mission planning and for real-time 
navigation, guidance, and control of the spacecraft for 
lunar missions. It is also of scientific interest since 
knowledge of the gravity can be correlated to physical mass 
concentrations. A number of gravity fields have been obtained 
from these data. However all have fairly poor prediction 
qualities. The object of this investigation is to develop 
a method which can be used for more accurate determination of 
lunar gravity parameters. 

DATA SOURCE 

Since 1966 the United States has successfully 
injected eleven satellites in close lunar orbits. The length 
of each of these parking orbits is given in Table I. 

Data was obtained using the three stations of the Deep Space 
Network (DSN) during the Lunar Orbiter flights, and by seventeen 
stations of the Manned Space Flight Network (MSFN) during 
the Apollo missions. The basic data types gathered by each 
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TABLE I 


SATELLITE 

PARKING ORBIT TIME 

LUNAR ORB I TER I 

76 DAYS 

LUNAR ORBITER II. 

335 DAYS 

LUNAR ORBITER III 

243 DAYS 

LUNAR ORBITER IV 

7 6 DAYS 

LUNAR ORBITER V 

179 DAYS 

APOLLO 8 

20 HOURS 

APOLLO 10 

62 HOURS 

APOLLO 11 

62 HOURS 

APOLLO 12 

86 HOURS 

APOLLO 14 

68 HOURS 

APOLLO 15 

144 HOURS 


tracking station consist of Doppler frequency shift and 
range data. Since using the range data requires an extremely 
accurate knowledge of the lunar ephemeris , this data type 
is not used for lunar orbit determination or selenodesy. The 
Doppler data is relatively insensitive to lunar ephemeris 
errors; hence it is the only data used. 

Two types of Doppler data are used: two-way and 

three-way data. For the case of two-way data a signal is 
transmitted by a station, frequency shifted by the satellite, 
and retransmitted to that same station. In the case of three- 
way data the signal is transmitted by one station (master 
station) , received and frequency shifted by the satellite. 


and retransmitted to the master staions and to any other stations 
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(slave stations) which are in the line-of-sight of the satellite. 

Essentially, the Doppler data type measures the relative velocity 

between the station and the satellite. No data are acquired 

when the satellite is occulted by the moon. The tracking 

stations use very sensitive double- conversion superheterodyne 

automatic-phase- tracking receivers. The transmitter, receiver, 

and station timing are driven by very stable atomic frequency 

standards. The stability of the Doppler frequency-shift data 

is a direct function of the frequency standard and is on the 

12 

order of one part in 10 

All the tracking data gathered to date covers a 
restricted region of the lunar surface. The inclinations of 
the Lunar Orbiters were limited to two regions, low inclinations 
between 10° and 21°, and a high inclination of about 85° (see 
Figure 1) . The inclinations of the Apollo missions ranged 
from about 178° to 154° (retrograde orbits) . (See Figure 2) . 
Hence there exists a large gap in the tracking coverage in the 
region between 26° and 85°. Any anomalous areas in the lunar 
mass distribution existing in the untracked region produce only 
secondary variations in the orbits that have been tracked and 
consequently are very difficult to determine accurately. 

Attitude control maneuvers performed by the satellite 
have damaging effects on the overall usefulness of these data 
for selenodesy purposes. In the case of Lunar Orbiter the 



LUNAR ORBITER IV 
(I = 85°) 



FIGURE 1 - LUNAR ORBITER ORBITAL GEOMETRY 
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attitude control system was not coupled in either the pitch 
or yaw direction. During all of the early Apollo missions 
(8, 10, and 11) a large number of attitude control maneuvers 
were performed using only one thruster quad. In each case 
the effect of these maneuvers was to produce not only the 
desired rotation but also a translational velocity for the 
satellite. The trajectory of the satellite is noticeably 
disturbed by these propulsive maneuvers. As pointed out by 
Lorell'*', in theory these effects can be modeled; but in 
practice it is a very difficult and costly operation. Further, 
because of engineering uncertainties and incompleteness of 
telemetry records, the reliability of results would be question- 
able. Every precaution possible is taken during the course 
of this analysis to use tracking data which is of free-flight 
quality (thrust free) . 

REVIEW OF PREVIOUS METHODS 

Two general methods have been used to obtain lunar 
gravity fields; the direct method and the indirect method 
or long-term selenodesy method. 

The direct method, principally used by Langley 
Research Center, attempts to estimate the gravity field 
from the Doppler tracking data. In this method the solution 
parameter set is augmented to include both dynamical state 
variables and lunar harmonic parameters. Large batches of 
tracking data covering long time periods from many satellites 
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2 3 

are used. Results obtained by Michael , Tolson and Gapcynski , 

4 

and Gapcynski, Blackshear, and Compton , are representative 
of this method. 

In the indirect method solutions are obtained using 
a two-step process. The first step is to regress Doppler ob- 
servations spanning about one day using an assumed gravity 
field to determine the satellite state at some particular time. 
The osculating states generated are averaged for each satellite 
period to obtain a mean value of the Kepler elements for that 
particular orbit. This process is repeated for each available 
day of tracking data. The second step in this method involves 
fitting the long-period perturbation equations to the Kepler 
element variations. The lunar gravity determination is 
accomplished by numerically integrating the perturbation 
equations and differentially correcting a set of gravity param- 
eters so as to give a best fit to the Kepler element variations. 

. 5 

The results obtained by Blackshear, Compton, and Schiess , 

6 7 

Lorell , and Risdal , demonstrate this method. 

OUTLINE OF THIS INVESTIGATION 

A new method for determining the lunar gravity 
field is developed which uses an empirical orbit determination 
(OD) and the long-period perturbation equations. The gravity 
determintion is a two step process as in the indirect method. 

The empirical OD. represents the Kepler orbital 
elements of the satellite as a six-dimensional time series. 
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The functional form chosen for each orbital element is that 
which best represents the long-period and secular variations 
in the orbital elements. In the first step of the process, 
tracking data is reduced to obtain simultaneous estimates 
for the orbital elements and the orbital element rates. This 
process is repeated for as many different spans of data and 
different satellites as are available (see Figure 3) . In the 
second step, the perturbation equations, linear in the gravity 
parameters, are solved using the orbital elements and element 
rates to obtain a field (see Figure 4) . 

This method differs significantly from the long- 
period approach that has been previously used in that the 
Doppler data reduction uses no assumed model and directly 
estimates a mean orbital element state. Further and perhaps 
most importantly, the simultaneous estimation of Kepler 
states and rates makes each solution independent of the next 
with respect to the gravity parameter estimator. 




FIGURE 3 - EMPIRICAL ORBIT DETERMINATION 





PERTURBATION 
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FIGURE 4 - DETERMINATION OF GRAVITY PARAMETERS 




CHAPTER II 


DYNAMICAL FORMULATION 


The motion of a close lunar satellite is governed 
by a perturbed Newtonian gravitational law of attraction. 

The perturbations arise from the non-central properties of 
the lunar mass distribution and from the disturbing effects 
of the gravity fields of the earth, the sun, and other planets. 
These perturbing attractions are small; their effects are at 
least one thousand times less than the inverse square attraction 
of the moon. 

The general equations of motion for a satellite 
referenced to an inertial moon centered coordinate system is: 

d 2 r . r , — /o t \ 

— 2 ■■ H — + a d (2 - 1) 

dt r 

where r is the position vector, is the product of the lunar 
mass and the universal gravitational constant, and a^ is the 
sum of all the perturbations. Although the perturbations are 
small numerically, their integrated effects are non-negligible 
and produce changes in the satellite's orbit. 

The purpose of this chapter is to present the 
basic equations of Newtonian or Keplerian motion. These 
equations will be extended to account for a primary gravity 
field of arbitrary mass distribution and to account for the 
perturbing effects of third bodies (e.g., the earth and sun). 
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NEWTONIAN MOTION 

The equations of motion for a satellite located 
at position r in a moon centered inertial system under the 
influence of an inverse square (Newtonian) gravity field are: 

r = ”^"3 ( 2 . 2 ) 

These equations possess closed form solutions which characterize 
particle motion as follows: 

1. The particle moves in a plane passing 
through the center of mass and has a 
constant angular momentum h = r x v 
which is perpendicular to this plane. 

2. The path of the motion is a conic section. 

Since the equations of motion (2.2) are a set of 

three second order differential equations, a complete solution 
is uniquely specified by six integration constants. Geo- 
metrically, these constants are usually interpreted as the 
classical Kepler elements given below (see Figure 5). 

a - semi-major axis of the orbit 
e - eccentricity of the orbit 
I - inclination of the orbital plane 
n - longitude of the ascending node 
ui - argument of pericenter 
T q - time of pericenter passage. 
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The variables a and e define the shape of the conic, 

I, £2 and a) are the Euler angles which specify the orientation of the 
orbit plane and axes, and T o locates the position of the satellite 
in the orbit relative to pericenter. The time of pericenter 
passage is relatable to the mean (M) , true (f) , and eccentric 
(E) anomalies using the following conic relationships: 


M(t) 



(2.3) 


M ( t) = E - e sin E (elliptical orbit) (2.4) 


f = 2 



(2.5) 


The satellite state is defined as a six vector of position 
and velocity; 


x = 


r 

v 


(r = v) 


( 2 . 6 ) 


Correspondingly the satellite Kepler element vector is a 
six vector of orbital elements 


a 


k = 


e 

I 

Q 


(2.7) 


u) 

M 
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The Kepler elements at any instant of time are non-linear 
functions of the satellite state and vice versa (see 
Appendix A) : 

x = x(k) (2 . 8 

k = k (x) 


It is important to note here that the entire Newtonian develop- 
ment is based on an inverse square force law which assumes a 
spherically symmetric gravitational potential of the form: 


U (r) 


r 


(2.9) 


Potential used here has the opposite sign to that used for 
potential energy in physics. 

PERTURBED MOTION 

If the satellite is now assumed to be under the 
gravitational influence of a primary planet (moon) of 
arbitrary non-homogeneous mass distribution, and of other 
planets, the gravitational potential of the satellite 
has the form: 

_ Va _ 

U(r) = + R(r) (2.10) 

where R is the disturbing function arising from the non- 
central effects. The equations of motion for the satellite 


now become: 



16 


r = vu = -y. -3 + R (2.11) 

^ r 9r 

g 

where a. = — R and V is the gradient operator. With R ^ 0 
9 r 

the Kepler orbital elements in general become functions of 
time. Since the non-central portion of the gravity field is 
small in comparison to the central part, the magnitude of 
variation in the Kepler elements is also small. The three 
second Order differential equations (2.11) describing the 
perturbed satellite motion can be transformed into the orbital 
coordinate system using the method of variation of arbitrary 

g 

constants to obtain the time rate of change of the Keplerian 

elements (perturbation equations) . The result is that the 

satellite state is now specified by six first order non-linear 

9 

perturbation equations : 


da 2 9 R 

dt na 9M 


de 

- V.l-e 2 Vi_ c 2 

9R 9 R 1 

dt 

2 V1 e 

9M 9 a) 


na e 

J 

da) 

_ Vl-e 2 [l 9R 

cot I 9 R 

dt 

2 e 9 e 

. 2 9 1 


na L. 

1-e 

df2 

csc I 9R 


dt 






( 2 . 12 ) 


(2.13) 


(2.14) 


(2.15) 
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dl 

dt 


dM 

dt 


esc 1 

27T 3 

na 'Vl-e 


cos 





(2.16) 


(2.17) 


where the 
equations 


mean motion 


"•V5- 


In vector notation these 


have the following general form: 


|| = f(k,t) (2.18) 

In order to carry out a solution for the equations, the 
disturbing function R must be specified and expressed in 
orbital plane coordinates. The disturbing function assumed 
in this investigation will have components arising from the 
non-central part of lunar gravity, R(^, earth and sun pertur- 
bations, R and R , and solar radiation pressure, R„^. 

Hence , 

R = R/. + R + R + R or) (2.19) 

© © SR 

LUNAR POTENTIAL FUNCTION 

The general form accepted by the International 
Astronomical Union 1 *^ in 1962 for the representation of the 
lunar gravity field is the spherical harmonic expansion 
expressed in selenographic coordinates. The selenographic 
coordinate system is moon fixed and oriented such that the 
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x-y plane is the lunar equatorial plane, the x-axis coincides 
with the mean earth-moon line through Sinus Medii, and the 
z-axis is directed along the moon's mean spin axis. The 
selenographic axis system rotates about the z-axis of the 
inerti’al system with a sidereal period of 27.32 days. 

The lunar gravitational potential is the solution 
to the Laplace equation (V U=0 ) in spherical coordinates and 
is expressed as follows : 


— R 1 


U (r , <f> , A ) = 


1 + Y. Z (t) 

1=2 m=0 v l 


% (sin<f>) < C &m cos m A + 


S. sin m A 
£m 


} 


( 2 . 20 ) 


In this equation R is the mean radius of the moon, P m 
is the associated Legendre polynominal, (r,<j>,A) are the 
selenographic radius, latitude, and longitude, and C , S £ m 
are the gravitational parameters which describe the non- 
central features of the lunar field. The objective of this 
investigation is to determine a finite set of (C„ , S„ } 

which accurately describes the lunar field. It is assumed 
that the origin of coordinates is located at the center of 
mass of the moon. Using this assumption, the terms of first 
degree (A=l) are omitted from the expansion since they repre- 
sent center of mass displacements in each dimension. The terms 
of degree A and order m=0 are known as zonals and are symmetric 
about the z-axis. Terms of degree A and order m are called 
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tesserals (£=m are sectorials) and are functions of all three 
dimensions. The potential of the moon is related to its 
density distribution D as follows : 


U(r) = 




dr 


( 2 . 21 ) 


where G is the gravitational constant and r' is the distance 
to the mass element. 

Two examples of the perturbations in the potential, 
one for an even degree zonal, C^q , and one for an odd zonal, 
C^q / are shown in Figure 6. 

The disturbing function for the lunar potential is: 


= U(r, <f> , A) 



( 2 . 22 ) 


where U here is given by equation (2.20). The disturbing 
function can also be represented as the sum of zonal and 
tesseral terms. 


R 


5 . 



(2.23) 


where each particular term, U^ m , is 
R 1 

p m 

u m “ p *. <sin + )<c « cos m 1 + s tm sin m x) (2 - 24) 



- 20 - 


Z 




EQUATOR 


PEAR SHAPE C 30 >0 


FIGURE 6 - PERTURBED POTENTIALS 
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A general formulation for transforming each term 
in the disturbing function of a primary body has been 
developed by Kaula^. During the course of this derivation 
the spherical harmonic potential is converted from the 
selenographic spherical coordinate system (see Eq. 2.20) 
to the Kepler elements referenced to the inertial 
system. The basic assumption in this derivation is that 
the inertial and selenographic systems share a common 
equatorial plane and that the selenographic system rotates about 
the polar (z) axis of the inertial system (see Figure 7) . 


Z 



FIGURE 7 - ROTATING COORDINATES 
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The disturbing potential expressed in orbit plane 
coordinates is : 


R*(k) 


oo 

10 - Z 


N-frr p (I,e)S (2.25) 

£=2 riFD 11 a* +1 


where 


& +! 


m /y 

P(I,e) - £ £ 


p=o q=-« 


Jimp Jlpq 


(2.26) 


and 


F « P (I) 


- Z 


(2l-2t ) i 


t t! (Jl-t) J (Jl-m-2t) ! 2 


( 2 A— 2 1 ) 


. £-m-2t _ 
sin I 


t (I) 


s=o 


cos s I £ 
c 


( Jl-m-2t+s\ / m-s \ 
c / \p-t-c/ 


1) 


c-k 


(2.27) 


(p, q, t, s, c, and k are dummy indices) 

where k is the integer part of • t is summed to the 

lesser of p or k, and c is summed over all values which 
make the binomial coefficient nonzero. 

The function G (e) in equation (2.26) is defined 
as follows : 


G (e) = (-1) 
Jlpq 


(i+e 2 ) V q l 


CO 

z 

k=o 


TO 6 

Jlpqk v Jlpqk 


2k 


(2.27) 


where 


e 


l+/l-e 2 


(2.28) 
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' &pqk 


£ eo ' 


h = k+q' , q'>0; h=k, q'<0 


and 


Q 


£pqk 


■ £ (;r) * 


h=k , q'>0; h=k-q', q'<0 

p'=p, q'=q for p_<£/2 

p'=£-p, q'=-q for p>£/2 


Also in equation (2.25) 


£.mpq 


'£m 


L S £m J 


£-m even 

cos [ U-2p) u>+ (£-2p+q)M + m(ft-e)] + 
£-m odd 


-i£-m even 


£m 


L Am J 


sin [ ( i,-2p) u> + (i,-2p+q)M + m(ft-0)] 
£— m odd 


(2.29) 


(2.30) 


(2.31) 


and 0 is the angular displacement between the inertial and the 
graphic axis systems. 
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The effects of the U^ m disturbing term can be 

. . 12 
approximated by using standard perturbation theory concepts 

Essentially each term gives rise to short-period and either 

long-period (and very long-period) or secular variations in 

the orbital elements. As applied to this dynamical system, 

short-period variations are those which are periodic with a 

period of an integer multiple less than the satellite orbit 

period. Long-period variations are again periodic with a 

period which is an integer multiple less than the period of 

the moon (28 days). Very long-period variations have a 

period less than the period of the peri focus rotation. 

Secular variations are nearly linear in time. Examining 

equation (2.31), these different types of variations are 

associated with the following factors : 


1. 

Short-Period : 

U-2p+q)M 

2. 

Long Period: 

m (ft- 0 ) 

3. 

Very Long-Period: («,-2p)w 

4. 

Secular: (£- 

2p) = U-2p+q) 


. 13 

Analysis has shown that the long-period satellite 
dynamics can be accurately determined from Doppler tracking 
data. This form of the perturbation equations , the long-period 
equations, can be obtained by averaging the lunar disturbing 
function, R , with respect to the mean anomaly. 
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R <r 27 / h™ 


(2.32) 


This calculation was also carried out by Kaula and the result 
simplifies (2.25) to the following: 


00 X, 

L t 


£ = 2 m=b ^ a^ + 1 P (I ' e) S ^mpq (to ,n , 0 ) 


(2.33) 


where q = 2p-Jl only and 


P’ (I ,e) = / F„ (I ) G (e) 

* m P apq 


(2.34) 


F £mp ( I ) is given by (2.27) and G ^ pg (e) is: 


G r, — ( e ) - 


( 1-e' 


1 ^r 1 ( 1-1 V 2a+e_2p i 

! ) d=o ^2d+£-2p'yy d J \ 2 ' 


in which 


(2.35) 


P' = P 


P' = *“P 


P 1 1/2 
P > a/2 


and S lmpq (t ^' 0) is 


s * 

2,mpq 


«,-m even 


cos [ ( £-2p) u) + m(fi-e)] + 


£-m odd 
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£-m even 

sin [ ( £ — 2p) to + m(fi-e)] 
£-m odd 


(2.36) 


It should be noted here that this entire long-period averaging 

14 2 

procedure is accurate to the first order . (Terms of order e 
and above are neglected. ) The long periodic variations experi- 
enced by R| can be either long-periodic in the ascending node 
(m ^ 0) , very long-periodic in the perifocus (£ - 2p ^ 0) , both, 
or secular (m = 0 , £ = 2p) . The secular variations are associated 
only with zonal perturbations of even degree xiJ . (Very loncr- 
periodic variations are also experienced by these terms.) Zonal 
terms of odd degree have only very long-periodic variations. 
Tesseral and sectorial terms have both periodic and long-periodic 
variations . 

Since the disturbing function, , is independent of 
the mean anomaly, two of the perturbation equations (2.12 and 
2.13) simplify to: 

= 0 (2.37) 


de 

dt 



(2.38) 


Therefore, to first order exactness, the semi-major axis of 
the orbit is a constant of the motion. 

The complete analytical equations specifying the 
long-period dynamics of a satellite under the influence of 
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a moon of arbitrary mass distribution are given by (2.37) , 
(2.38), (2.14), (2.15), (2.16), (2.17), (2.33), (2.34), (2.27), 

(2.35), and (2.36). As an example of the form of long-period 


equations, the perturbation equations for the C2Q zonal term 
and the C ^ sectorial are presented: 


For C 


20 


da de dl^ n 

dt dt dt 


(2.39) 


0 n C or . /R \ 
dfi _ 3 20 I _e ] 

dt 2 , . 2 v 2 \ a I 

(1-e ) \ ' 


cos I 


(2.40) 


doj _ 3^ 
dt 4 


{-) 

(l-e 2 ) 2 V a / 


(1-5 cos I) 


(2.41) 


dM 

dt 


n f ‘ t 7T- 


'20 


2, 3/2 
e ) 


(3 cos^I-1) 


(2.42) 


For C 22 : 


da _ de 
dt dt 


(2.43) 


dl 

dt 


3n 

(l-e 


'22 


2 ^ 2 



sin I sin2(fi-0) 


d fl 

dt 


3n 

(l-e 


'22 


2,2 



I cos2(fi-0) 


(2.44) 


(2.45 
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dw 3_ 
dt 2 


dM _ 
dt “ 


n 


'22 


, , 2.2 
(1-e ) 




2 


2 


cos2 (n- 0) (5sin I - 


• 2 _ 
sin I 


(1-e 2 ) 


cos2 (ft-e) 

372 


— 


2 ) 


(2.46) 


(2.47) 


For both sets of perturbations the dynamical equations are 
coupled non-linear equations (linear in the gravity parameters). 
A principal difference between the zonal term and the sectorial 
(and this is true in general) is that zonal perturbation 
equations are autonomous whereas the others are not. 

EARTH AND SUN PERTURBATIONS 

The disturbing function for a third body (earth or 
sun in this case) located at r^ on a satellite located at r 
has been computed by Brouwer and Clemence^^. 


R 


3 




3, 

COS ip 


— 

j + ... 


+ 

(2.48) 


where represents the third body mass and ip is the 
angle between r and r^ (see Figure 8) and r^ >> r. This 
expansion is accurate to third order. In order to 
obtain the long-period portion of this disturbing function, 
it must be averaged over the satellite orbit. 
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The averaging with respect to the anomaly of this 

1 7 

disturbing function has been performed by Liu and Lorell 
and has the following form: 


where 



(2.49) 


F = [JlV I 

2 rj 2 


- (1+3/2 e 2 ) + 3A 2 (2e 2 +l/2) + (1-e 2 ) 


i 


(2.50) 


F 3 = 


Ae [ 3 


Ae 13(1+3/4 e 2 ) - 5A 2 (e 2 +3/4) - ^ B 2 (1-e 2 ) 


(2.51) 


] 


and A=u'P,B=u*Q. u is the unit vector in the direction 
of the third body, P is a unit vector to the perifocus and Q 
is a unit vector in the orbit plane orthogonal to P. 

Hence, the total dynamical effect on a satellite 
in lunar orbit due to either the earth or sun is formed by 
evaluating the perturbation equations using (2.49) as the 
disturbing function. 

PHYSICAL LIBRATIONS OF THE MOON 

In the derivation of the lunar disturbing function 
it is assumed that the selenographic axis system to which 
the gravity parameters are referenced rotates about the polar 
axis of the inertial system. In addition to this polar 
axis rotation, the selenographic system undergoes additional 
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small rotations about all three axes due to precession and 

physical librations. Since these rotations are not included 

in the analytical formulation a correction must be applied 

18 

to the long-period equations. Analysis shows these effects 
can be eliminated if the selenographic system fixed at epoch 
is used. 

SOLAR RADIATION PRESSURE 

The perturbing effects of solar radiation pressure 

19 

have been derived by Kozai . The variations in the orbital 
elements are : 

Aa SR = 0 (2.52) 

(first order approximation) 


Ae gR = - F A*Vl-e 2 cos 2E+B*(^ sin 2E-2esin E+|e ) 

^ an 


E, 


(2.53) 


FP* Z G. 3 

AH' = r~ ■- . .. cos 0 ) [(1+e ) sin E - -r sin 2E - ^ eE] + 

SR a nVl-e 2 4 2 


(2.54) 


2 - - r ^ e 


e sin a) [cos E - ^ cos 2E] 


Aft 


FC* i 2 e 

er> n 7— — T I sina) [(1+e ) sin E - 7 sin 2E - ^eEl - 

SR — sinfvi-e 2 4 2 

: (2.55) 


an 


Vl-e 2 


cosui [cosE - -j cos2E] 
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A co 


-F 


SR a r 2 

a n e 


A*Vl-e 2 (e sin E + 7 sin 2 E - ^ E) + 


1 
4 


3 

2 


E. 


(2.56) 


B* (e cos E - 2 cos 2E) 


- COS I Aft 


SR 


E. 


AM SR = I A * [ ( 1+e 2 ) sin E - 

a n 

B* Vl-e 2 [cos E - | cos 2 E] 

where F is the magnitude of solar acceleration and A* = P • u @ , 

B* = Q • u , C* = W • u . u is a unit vector in the direction 
® © © 

of the sun, P, Q, W are vectors to peri focus, in the orbit 
plane orthogonal to P, and normal to the orbit plane orthogonal 
to P and Q, respectively. 


| sin 2E - |eE] - 

fa 


- V 1 -e 2 ( A 8 __ + 


( A ftg R + cos I Aft SR ) 


(2.57) 


F 



(2.58) 


<\j 

where A is the satellite effective cross-sectional area, m is 

the mass, S is the solar constant, r the distance to the sun, 
© © 

y. is the reflection coefficient, and c is the speed of light. 

These perturbations are only acting on the satellite 
when it is in sunlight. A test used to determine the sun's 
visibility is (see Figure 9) : 




34 


where . 


and 



i lei 
> | 6 | 


perturbations are zero 
perturbations are present 


B| = |tan -p| 

(2.59) 

a | = Icos' 1 

(2.60) 


1 = r - r 
© 

(2.61) 


SUMMARY 

The perturbation equations for the long-period 
motion of a lunar satellite under the gravitational influence 
of a primary of arbitrary mass distribution and the earth, 
sun, and solar radiation are presented. The general form 
of these dynamical equations is as follows: 


da 

dt 


0 


(2.62) 


de 

dt 


dl 

dt 


E 

£ ,m 


de» de^ de de er . 

t,m . © . ® . SR 

r -i i ■ -i i ' 


dt 


dt dt dt 


£ 

£ ,m 


dl. dl dl dI C r, 

£ ,m , © , ® , SR 

— q-— ~T — 


dt 


dt dt dt 


E 

£ ,m 


£ ,m © ® SR 

dt dt dt 


(2.63) 


(2.64) 


do 

dt 


dt 


dt 


dt 


(2.65) 
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doo 

dt 


dM 

dt 


Z 

i ,m 

Z 


l ,m 


£ ,m 


© 


dt 


dt 


+ 




SR 


dt 


dt 


dM 


£ , m 


dt 



dM . , dM SR 
dt dt 


( 2 . 66 ) 


( 2 . 67 ) 


These long-period perturbation equations form an 
analytical basis for the development of the empirical OD. 



CHAPTER III 


EMPIRICAL ORBIT DETERMINATION 


One method of modeling the perturbed motion of a 
lunar satellite is to represent each of the Kepler elements 
of the orbit as an independent time function . This repre- 
sentation forms the basis of the empirical orbit determination 
method. Using Doppler tracking data and a minimizing process, 
estimates are obtained for the parameters which characterize 
these time functions. 

The purpose of this chapter is to present the 
theory and equations used to represent the long-period motion 
of the lunar satellite. 

MATHEMATICAL THEORY 

The long-period time dependence induced in each of 
the Kepler elements by non- central gravity perturbations is 
given by the perturbation equations (2.62-2.67). 


dk 

dt 


f (k,t) 


Since these equations are non-linear, general closed form 
solutions are not obtainable. The non-central effects are 
extremely small compared to the central body term; consequently 
solutions can be approximated. 
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If an analytic quadrature is performed on each 
perturbation equation, a set of integral equations results: 

k(t) = k(t ) +f f (k , t) dt (3.1) 

° K 

The kernels, or forcing functions, appearing above are non- 
separable, non-linear functions. If only the perturbations 
arising from the lunar disturbing function are considered, 
these kernels can be categorized into two types : 


1 . 

f (k ,t) 

= f (k) (autonomous) 

(3.2) 

2. 

f (k ,t) 

= g(k)sin myt + h(k)cos myt 

(3.3) 


m = 1 , 2 , 3 ,... (perturbation order index) 
where y (1/27.32 rev/day) is the rotational rate of the moon 

(0 = Yt) . (This categorization does not apply when the 
standard Kepler elements become singular, for nearly circular, 
nearly equatorial orbits. However, a similar approximation 
can then be used with a modified orbital element set. ) The 
first or autonomous kernel corresponds to zonal perturbations, 
and the second to tesseral and sectorial perturbations. If it is 
assumed that for periods of time of 24 hours or less the 
magnitude of variation in the Kepler elements is small, then 
f(k), g(k) , and h(k) can be considered constant. 

Solutions to equations possessing autonomous kernels 
lave the following simple form: 
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k(t) = k(t o ) + f(t o )(t-t Q ) (3.4) 

These solutions have the linear properties of secular variations. 
Using the same time approximation made previously (as for 3.4), 
solutions for the non-autonomous kernels can be given as follows: 

k(t) = c n (t )+ c„ (t )sin m y(t-t ) + c_ (t ) cos m Y(t-t ) 
lo2o o3o o 

(3.5) 

where c , c~ , and c, are constants and c (t 1 = kCt ) ^ cT Ct ) , 

-L ^ -L O O JO 

If the time period of the data span (t-t Q ) is 
smaller than the period of a spherical harmonic perturbing term, 
this solution can be expanded in a truncated Taylor series : 

k ( t) = c‘ 1 (t Q )+ c 2 (t Q )£m y(t-t Q ) + ... 



or 

k ( t) = I + A.t + l„t 2 + 

o 1 2 

where the Aj represent vector constants. 

Hence the functional form for the Kepler elements 
which best represents the long period and secular effects 
is that given by (3.7). Two typical Kepler elements (for 
example ft(t) and e(t)) are represented as follows: 


•] 

(3.6) 

! 

(3.7) 
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ft(t) = ft + Q,t (secular variation) 

o 1 

2 

e(t) = e + e^t + e^t (secular and long- 
° period variation) 

The terms fl , e Q , e^, and e ^ are examples of Keplerian 

parameters determined by the empirical OD method. 

Since the third body and solar radiation perturbations 
affect the satellite state in orbit, these effects must be 
included in the Kepler element equations. The actual six- 
dimensional time series used to represent the Keplerian 
motion of the satellite is: 

a(t) = a constant (3.8) 

e ( t) = e + e.t + e 0 t 2 + 6e + 6e + (3.9) 

o 1 2 © © SR 

I(t) = I + I.t + I„t 2 +61 +61 + 61 (3.10) 

O 1 A © © bR 

J2(t) = fi + fl.t + ft n t 2 + 6fi + 6fi + 6^-,^ (3.11) 

o 1 2 © © SR 

2 

a) ( t) =co + u),t + u) n t 4* 5 a) + 6 a) + 6a )«_ (3.12) 

O 1 2 © © S R 

= M + M,t + M„t 2 
o 1 2 


M(t) 


(3.13) 
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where the orbital element variations 6k are found by numerical 
integration. The third body and solar radiation perturbations 
are those developed in the previous chapter. No explicit 
third body and radiation pressure perturbations are modeled 
for the mean anomaly. Hence the time-series for M(t) includes 
the perturbations of the moon, earth, sun, and solar radiation. 
The reasons for this representation are associated with the 
semi-major axis estimation and are discussed later in this 
chapter. 

Terms of quadratic order are used in the time-series 
since for the time periods (24 hours or less) over which 
solutions are obtained, these will adequately (by conservative 
standards) represent long-period perturbations of up to order 
seven (m=7) . (A harmonic perturbation of order seven has a 
period of four days. Any one day segment of this harmonic can 
be well represented by a quadratic function.) 

DATA REDUCTION 

# 

The Doppler observation (p) is a scalar quantity 

that is a non-linear function of the satellite state, k, 

the relative earth/moon configuration and the earth-based 

tracking station position and rotational velocity. During 

the orbit determination only the estimates of the satellite 

state are refined. The tracking data P obg is related to the 

21 

satellite state as follows (see Appendix B) 
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P 0 bs = + n (3.14) 

where n is the random noise associated with the physical 
measurements. The measurement errors are assumed to have the 
following properties : 

Exptn] =0 (3.15) 

Exp[n 2 ] = a 2 (3.16) 

2 

where Exp is the expectation operator and a is the variance 
of the Doppler noise. Since the Doppler is a non-linear 
function of the Kepler state, the output equation (3.14) must 
be linearized about a reference trajectory. 

If k(t) is the true Kepler state, k*(t) is the 
reference state, and Ak the deviation, then 


Using (3.14) 


k ( t) = k* ( t) + Ak (t) 


p ob s « + m Ak + 


(3.17) 


+ n (3.18) 


Defining 


and 


Ap (t) = P obs ~ P [k* (t) ] 

J(t) =■(! Si)* 
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the linearized relationship between the deviation in the 
Kepler state and the deviation in the Doppler is : 

Ap (t) = J Ak (t) + n (3.19) 


Using a similar procedure, a linear relationship 
relating the Kepler state k and the parameter vector K can 
be obtained. K is an (n*l) vector consisting of the solution 
parameters 


K = 


a 


e 

o 



I 

o 


(3.20) 


M, 


Since these Keplerian parameters are constant over 

the trajectory, an expansion of K can be performed about a 

reference set K* at some time t . 

o 


Since 


K(t Q ) = K* (t Q ) + AK(t Q ) 
k = k[K(t Q ) ,t] 


(3.21) 

(3.22) 


k ( t) = k[K*(t Q ) ,t] 



AK ( t ) + 

o 


then 



43 


or 



(3.23) 


Deviations in the parameter set AK(t Q ) are then linearly 
related to deviations in the Kepler elements. The output 
Doppler equation (3.19) can now be expressed in terms of the 
solution parameters. 


Ap(t k ) = H(t k ) AK(t Q ) + n (3.24) 


where 



The primary objective in orbit determination can be 
simply stated as follows: Given a batch of Doppler measure- 
ments valid at many times t^ , , . .., t r , find the set of 

Keplerian parameters which best fits these data. The linear 
relationship between the Keplerian parameters and the Doppler 
can be generalized as follows : 

Ap = [H]AK + TT + s (3.25) 

where 

Ap is the (r x l) observation vector 
H is the linearized set of functions 
relating the observation and state 
(r x n) matrix 

AK is the column vector of Keplerian 
parameters (n x 1) 
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n is an (r x l) vector of observation noise 
s is an (r x l) vector of systematic modeling 
errors . 

Since the Keplerian parameters K only attempt to 
represent the long-period dynamics of the orbit, then at any 
point in time there will always be a systematic residual (s) 
associated with the unmodeled short-period satellite variations 
Various data smoothing techniques were used on the Doppler data 
in an attempt to eliminate these short-period variations. None 
of these approaches has succeeded. Attempts to model the short 
period variations directly in the empirical OD also have not 
been successful. 

Given the linear system of equations (3.24) the next 

step in the OD process is to formulate an estimation scheme 

which minimizes the estimation error and yields a best estimate 

of the Keplerian parameters. If the random quantities in the 

dynamical system are assumed to be normally distributed, there 

22 

are three linear estimators (least squares, minimum variance, 
and maximum likelihood) which all would yield the same re- 
sults and could be used for this function. However, since 
the errors in the dynamical system are not random, and since 
the data to be processed is in batches , the weighted least 
squares estimator was chosen to perform the data reduction. 
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An error function, e, for r observations is given 
as follows: 

e = [ v^ v] (v=n + s) (3.26) 

or 

e = (Ap - HAK) T W ( A p - HAK) (3.27) 

where for one data type W is an (r x r) diagonal weighting 
matrix. For a minimum, the first variation of e must 
vanish ; 

<5e = 0 = - (Ap - HAK') T W H6K - (H6K)' I W(Aj - HAK) (3.28) 

where AK is the value of K that extremizes e and 6K is the 
variation in K. Since this is a scalar equation, 

0 = - (Ap - HAK) T W H 6K (3.29) 

But 6K is arbitrary, hence 

~ - T 

(Ap - HAK) W H = 0 (3.30) 

T 

Transposing and solving (W = W ) : 


1- T —1 T 

AK = [H W H] H W Ap 


(n x 1) 


(n x n) 


(n x 1) 


(3.31) 


(for minimum = 2 SK^ 1 [H^WH] ^ 6K > 0 for arbitrary 6K) . 

Since this minimization process was obtained by 
linearizing a non-linear set of equations, the least squares 
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estimation must be performed iteratively. The processing of 
r Doppler measurements and the resulting set of differential 
corrections AK constitutes one computing iteration. The 
convergence criterion for two successive iterations is as 
follows : 


£ Ap(t.) 2 
i=l 1 

x a p(V 2 

■ 1=1 


(K-l) 


(K) 


1 < 6 


(3.32) 


-4 

where 6 = 10 and (K-l) and K designate computing iterations. 

The entries in the H matrix are the partial derivatives 
of the Doppler with respect to the Keplerian parameters. These 
terms are found by using the chain rule for differentiable 
functions : 

= ii • (3.33) 

9K 9x 9k 9 K 


Each of the derivatives used in the above equation is given 
in Appendix C. 

SEMI-MAJOR AXIS DETERMINATION 

Studies using both pseudo and real Doppler data 

have shown that the least squares process does not converge 

when the semi-major axis is included as an independent parameter. 

This condition was also observed by the Smithsonian Astro- 

23 

physical Observatory in earth satellite work 
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Since the mean motion of the orbit, even in the 
presence of perturbations , is nearly inversely proportional 
to the semi-major axis to the three halves power, then the 
estimated mean motion can be used to imply a semi-major axis. 

The constraint equation between the mean motion and semi-major 

axis must include the long-period effects of the earth, sun, 
and solar radiation and some representation for the lunar field. 

The third body perturbations on the mean anomaly have been 

24 

given in the previous chapter. The LI lunar gravity model 
(see Table II) developed by Langley Research Center is used. 


TABLE II. LI GRAVITY FIELD 


20 

-2.07108 

X 

O 

1 

it* 

22 

0.20716 

X 

10~ 4 

30 

0.21 

X 

10 -4 

31 

0.34 

X 

i 

o 

i— i 

33 

0.02583 

X 

io -4 


The second degree terms in this model were deter- 
mined on the basis of astronomical observations of the moon's 
physical librations and hence can be assumed known to at 
least one signigicant digit. The third degree terms were 
obtained from previous analyses of Lunar Orbiter and Apollo 


data. 
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The average value of the semi-major axis (a Q ) with 
respect to mean anomaly is found by setting the mean anomaly 
perturbation equation equal to the empirically determined 
mean motion. Using equations (2.62) and (3.13), the following 
relationship is made: 


M = M 1 + 2M 2 t 


(3.34) 


M 




c„„ + M c„, + M c,„ + M c„ + M C„ + + M o + m sr 

20 22 30 31 33 


where and M 2 are the estimated Keplerian parameters. Since 
this is a non-linear equation in a o , it is solved using 
Newton's Method. The value a Q is the geometric average of the 
semi-major axis under the perturbing influences of the LI 
lunar gravity field, the earth and sun, and solar radiation. 

For the case of secular perturbations on an earth 
25 

satellite, Kozai has shown that a mean value of the semi- 
major axis a yields the average satellite position in the 
orbit. This value is derived such that the deviations in the 
position of the orbit due to perturbations averaged over the 
orbit yield only short-period variations. 



6r dM 


£ cos (c M + d ? ) 


(3.35) 


where 6r = r TRUE - r L0 ^ G p ERI0D and b ? and d^ are constants. 
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Hence the value a yields the proper mean position 
over the orbit. Kozai has derived the relationship between 
the average value, a Q , and the mean value, a, for the 
perturbation : 


a 




(1-3/2 sin 2 I) 
(1-e 2 ) 3 / 2 J 


(3.36) 


also since: 


M 


20 



C 2Q (1-3/2 sin 2 I) 
(1-e 2 ) 3 / 2 . 


(3.37) 


then the following modified Kepler law is valid: 


where : 
JJ 



y 



/ R e\ (1-3/2 sin 2 I) 

W (1-e 2 ) 3 / 2 j 


(3.38) 


(3.39) 


Hence the presence of the C 2Q term essentially changes the 
effective mass of the moon. 

Odd degree zonal, tesseral, and sectorial terms in the 

semi-major axis constraint equation do not require mean value 
2 6 

corrections . Analysis has shown that the correction for the 
sun and solar radiation are on the order of a foot or less; hence 
no factor is included for these perturbations. The mean value 
correction arising from the earth perturbations is included. 
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An analytical procedure was used to derive the rela- 
tionship between the mean and average values of the semi-major 
axis for earth effects. The basic formulation for the dis- 
turbing effects of the sun on the moon's orbit was used as a 
model . The secular effects on the mean motion generated by 
earth perturbations is as follows 


M 

© 



(3.40) 


where terms of order e are neglected and n is the mean 
motion of the earth. Extending the mean value solution 
given by Brouwer and Clemence to allow for varying inclinations 
(see Appendix D) , 


a 


a 

o 


n 9 

1 + -^2 (1-3/2 sin I) 


(3.41) 


l n -J 

2 

where again terms of order e are neglected. Both of these 
formulas agree well with those given for the sun's pertur- 
bations on the moon (case in which I = 23°). 

The general procedure for implying a mean value 
semi-major axis has two steps. First, the average value of 
the semi-major axis, a , is determined once per iteration 
using equation (3.34). Then, the mean value, a, is calculated 
using equations (3.36) and (3.41); 


t 


a * a o [1 + e 20 + 


(3.42) 


where the e terms are the mean value corrections. 



51 


The quantity a Q is the average value of the 
osculating semi-major axis. The quantity a is introduced 
to make long-period perturbation theory represent the average 
satellite orbit. Since this variable is introduced to insure 
compatibility between the long-period and associated rect- 
angular equations of motion, it is only used for data reduction. 

After convergence has been reached, the mean anomaly 
rate, M, is adjusted to remove the earth, sun, and solar 
radiation effects. 

= M. + 2M~t - M - M - JVL, (3.43) 

'i. 1 2 © © SR 

Here is the anomaly rate arising only from lunar gravity. 

The average semi-major axis value, a ^ , for the lunar gravity 
is found by solving the following equation: 



It is this value of the semi-major axis that completes the 
K parameter set and is used for gravity field determination. 

One such value of a^ is found for each solution. 

SOLUTION PARAMETERS 

The output from processing a batch of Doppler measure- 
ments is a best estimate for a set of Keplerian parameters, K. 
Since the third body perturbations are modeled separately in 
the OD process, and the mean anomaly parameter is adjusted 



52 


for third body effects, these Keplerian parameters represent 
the variation in the Kepler orbital elements arising from 
the non-central part of lunar gravity. These solution 
parameters give a simultaneous time history of the Kepler 
elements and element rates valid over the data span. A 
detailed block diagram of the empirical processor is shown 
in Figure 10. 

The time histories of the Kepler elements and 
element rates are used as input data to a second processor 
which fits lunar gravity harmonics to the Kepler element 
rates. Since the solution parameters provide continuous 
time functions for the orbital elements and rates they can 
be sampled at any desired frequency. The long-period lunar 
gravity effects have periods which are much greater than a 
typical lunar orbiter period (lunar effects have periods of 
days whereas a typical satellite period is three hours) , 
hence there will be no aliasing of gravity information if 
the element states and rates are evaluated once per satellite 
period. For example, if a typical data span contains five 
satellite periods, then five sets of Kepler elements and 
element rates are obtained. 

The Kepler element rates , which form the actual data 
for the gravity estimation, consist of the following five 


rates 



EARTH AND 
SUN EFFECTS 
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FIGURE 10 - ORBIT DETERMINATION BLOCK DIAGRAM 
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k 



(3.45) 


U) 

M 


The semi-major axis rate is zero in long-period theory, so 
it is not used. 

The five orbital element rates are simultaneously 
processed to obtain gravity coefficients. Since the input 
data consists of five different quantities (e , I , to , ft , M ) , 
a weighting matrix is required to define the relative accuracies 
of each of the rates. 

If the empirical OD method could model the Doppler 

such that the residuals remaining after convergence were 

T -1 

normally distributed random errors, the [H W H] matrix 

given by equation (3.31) would be the covariance matrix 

associated with the K solution set. Since only the long-period 

satellite dynamics are represented in the empirical process, 

T - 1 

then the [H W H] matrix is not the covariance matrix of the 

28 T -1 

process . However since the terms in the [H W H] matrix 

do reveal the relative sensitivity and correlations among the 

solution parameters, it is assumed for weighting purposes that 

these terms can be regarded as variances and covariances in 


the conventional manner. 
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The weighting matrix for the orbital element rates , 

29 

A, is a (5x5) matrix having the following form : 



2 

a • t . » 

e 12 

t 13 

t 14 

T 15 


2 

t 21 CT i 

• * • 

• • • 

• • 
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2 
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• 

• • • 

_ t 51 

• « • 

2 

o • 

U) 

• • • 

• 

• 

• 

2 

°M 


2 

where o£. are the error variances among the rates and t . . are 

K 1] 

the covariances. It is assumed that the mean error in the 
Keplerian parameters is zero. Each of the orbital element 
rates has the following form. 


& = n x + 2^ 2 t (3.47) 

The variance of the error in £2 is found as follows : 

a* 2 = Exp[(fi 1 - ^ 1 ) 2 + 4t(fi 1 - £2^ (^ 2 - £2 ) 

+ 4 1 2 (£ 2 2 - £ 2 2 ) 2 ] ( 3 . 48 ) 

or 

°5 2 = a l 1 + 4t oov,E n 1 i: a2 ) + 4t2 0 a 2 <3.49) 

= (SI 1 - op and = (0 2 - 0 2 ) . 


where e 
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The covariance terms among the rates (e.g., cov (e^e^)) 
are formulated in a similar way (see 3.48) 

cov (e*e*) = cov ( e e ) + 2t[cov(e e ) + cov(e„ e )] 

Q a) “]_ w 2 U 2 W 1 

+ 4t^ cov(e 0 c ) (3.50) 

a 2 “2 

Then if it is assumed for weighting purposes that 
T - 1 

the [H W H] is the covariance matrix, each of the entries 
in the A matrix can be found from the appropriate slots in 

T - 1 . . 

the [H w H] matrix. Hence a weighting matrix is auto- 
matically obtained for each set of orbital element rates. 

PSEUDO DATA SIMULATIONS 

In order to demonstrate the operational capabilities 
of the empirical orbit determination method, pseudo Doppler data 
were generated from numerically integrated trajectories and 
converged solutions were obtained. The estimated Keplerian para- 
meters from two such typical convergences are presented. In both 
cases a triaxial lunar gravity field is assumed (C 2 Q and C 22 
harmonics) and the earth, sun, and solar radiation perturbations 
are also included. No noise or biases are added to the pseudo 
data. 

The first data simulation was generated for a Lunar 
Orbiter V (polar) orbit. The data span contains tracking data 
from three stations (Goldstone, Calf., Madrid, Spain, and 
Woomera , Australia) and is approximately 21 hours 30 minutes in 
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duration. The epoch data and initial conditions for this 
orbit are : 

Epoch Date : Aug. 9, 1967 7 hours 20 min. 

Initial Conditions : (Selenographic) 

a = 8,324,332 ft. e = .27618984 

I = 84?764923 Q = 70?2050009 

to = 1J8616071 M = 244°73644 

The Doppler residuals (Ap) associated with this convergence are 
shown in Figure 11. These residuals are systematic and have the 
general form of the unmodeled short-period orbital variations. 
The residuals possess a mean of .0216 feet per second (fps) and 

a standard deviation of a = .1 fps. 

The Keplerian parameter set used consists of thirteen 

terms. The best estimates for these parameters are: 

e = .27600103 
o 

e 1 = .46568305 * 10 -9 
I Q = 84 ? 765987 

1^ = .89524714 x 10 ^ deg/sec. 

70 ?202817 

n = -.61366823 x 10~ 6 deg/sec. 
a 2 = .74160675 x 10 12 deg/sec 2 



-.12242802 x io 




oo o So 8| 


(Sdd) swnaisau uaiddoa 
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M = 244 ?73515 
o 

= .031390311 deg/sec. 

M 2 = -.55419592 x 10 -11 deg/sec 2 
The constrained values of the semi-major axis obtained are: 
a Q = 8,323,991 ft. 
a = 8,324,561 ft. 

In order to determine the quality of this converged solution, 
comparisons are made between the numerically integrated source 
trajectory and the solution obtained. Figures 12-15 show the 
variations in the six Kepler elements for both the converged 
solution and the source trajectory. The variations presented 
for the eccentricity (e) and the three Euler angles (I, Q, to) 
(see Figures 12 and 13) show the actual variations of these 
elements plotted on common axes. Since the variations in the 
mean anomaly are very large and the semi -major axis has a large 
magnitude, the differences in the converged solution and source 
trajectory values are shown for these variables (see Figure 14). 
As can be seen in Figures 12 and 13, the estimates of the 
eccentricity, ascending node, and perifocus have slight biases 
at epoch. The inclination parameter has a slight rate error. 

The mean anomaly and semi-major axis variations are only short- 
periodic and display essentially no growth characteristics. A 
plot of the difference in position between the source trajectory 
and the converged solution is shown in Figure 15. The 400 ft. 
bias between positions is relatable to the error bias in the 
estimated eccentricity parameter at epoch (e ) . The slight 
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FIGURE 13 - LUNAR ORBITER V ORBIT 
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DIFFERENCE IN POSITION AS A FUNCTION OF TIME 


FIGURE 15 - LUNAR ORBITER V ORBIT 
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time variations in the position differences arise from small 
errors in the eccentricity and mean anomaly rates. 

Analysis of the functional properties of the partial 
derivatives used in the least-squares convergence shows that 
two of the state variables have similar sensitivities to the 

data ; 


H * m (3 * 51) 

This lack of separability arises primarily from the poor 

geometrical configuration from which the Doppler measurements 

are taken. Plots of the functional behavior of the five basic 

partial derivatives used in the solution are shown in Figures 16-18. 

The presence of these nearly equal sensitivities among the state 

T -1 

parameters leads to high correlations in the (H W H) matrix. 

These correlations lead to linear combinations among the state 
parameters being estimated. For this particular solution very 
high correlations (p > |.9 | ) exist between the following sets of 
parameters : 

Correlation Coefficient Parameter Pair 


-.98 

I 

o o 

.97 

H 

0 

o e 

.91 

I l'“l 

-.97 

SO 

0 

■> 

e 

0 

-.92 

/ U) 


.97 



^ (FT/RAD-SEC) dp (FT/SEC) 


— 65 — 


!!!!!£!!!!!£ a ! BBaBBBB>aaaBBBaaaaBaaaaBBaaaaaaaBaaaBaaaaaaaaaaaaaaaaaaBB 

*■■■»■■»> »■■■■■■■■■> ^bbbb bbbbb bbbbb bbbbb bbbbb bi 

«■■■■■■■■■■ BBBB BBBBB ■■■■■■■ BBB6B Bl 


■■■•■■■■I 


■■■■■■■■■I 


mmmmmmmmmmmmmmmmmmmmmmmmmrnmummm mmuwxmmmummnmmi *bbbb ■■■■■ ■■■■■■■k# a a bi 

■ ■ a » ■ ■ i 

lamkiai 


■ a laniniiii iiiiiiaiiik 


ar laaaiaiiMimiaiiiiiaiapiiB 


BaBBBBBiK^iBBBaaaBarniBBBBBaBBaB^Bai 


■ IBBB jmBIBIBBBBBIBBBBBBBBI 


£?*•!££ !!!!!£!^S£!5 aaa 5 aL, * aaaaaaaaaaa-aaaaaaaaaariBaaaaa|laBIJ * aB 

■LaaflaBBBBBBB-jraBBaiBBB Bar v i ■■ a a * a ■ a * ? a ■ 


abb bbbbb b*bbbbbbb 


■maaaa' 


1BI 


’lBaraaBBaar . bb 


IT i -BB BABBBBBP A 1BBBP1BB 


BBBBB i UBBIIBI 
|BBBr:fl*]Bai£BBB| 


•■■■>■■■■» i iaa £ f «Wl IBBBrflBBBr fl' BBBBaBBBVillBBJtBBBB'Bt'BBBBJB 


iaaaa ■■■■• Ji !£Sif aaaaa ■•■■■■■■■■» ia aaaaBBBBB * bkbbbbbbbbv: ib^b bbbbb bb*ibb~ bbbbb ■ 


!££■■■» aa - .<BBaL!BBBB 'JflK.BBRBOBBBKIB.IBaarBBBB'lflMB^BiB 


[*BB*4BBfll£BBBBr4BfcBBBB71BI 


-1000 


-2000 


-3000 


■J !!!!! 1 !!!£!?■■■■■•■■■*■ ,B-BBajBBa,BBBajBB - aBBB »”u BB * BBB """i*Ml 

£2i!r , !2.;222?!2S225rJ222^222?!S v 2! BI,BPIBBBBa:4BBB:,BBBB,,BI,BBBfiBf:BBaBB ' ,BBBBBB ^fl 
b222i.! 7 iS!251?!!22f.!!!!i;!i!^!5 , !S B E ,B,, ' JB, * BBB,ai,:iBk:BB:B,BBBBBi:iBaaB,BaBBBB i| 

HllLiwAAMMAM£^MIMAMMMMBMMaa£aaIaaaaaaBaaiiaaa»«!!!*2''i!i5a!"ir!i 


■*a^a&' ar-aaaaa^iaaa^BvBaBBwaaaaiv. flBaaa'iaaa^B^aBBBPBBBaBajaHHHMSiMMMHHIM 
| BBBa i j BBBBBBBBB !- B il BBB ^|?!l a > : lll!g^^«»t?«-iB«BB >i«BBnnBiiiiDiSaS : .iSiiigiBirSB 
IUiM||||MMflMIBW^^^^^^^^M' :B » BB ' BBBBB »»^* 4 »»»»>«iH 


HUM 


■■?■■■» Jlllllllll ilBBBBOBBII 


i 2 s 5 : 5 s:ss:"-: 2 s 5 s 2 ! 55 t 2 £ss:s 252 s 5 "siiasisssssisiii-^asiasizssi_assssiiissi: 

If gBggggggggygBBBBBBBBBBlBBBBBBBBBrrf] 


BkdlBBBBBBBBBPiBBaaBBBIBBaaBBIjBBBBBMBBBi 


BBBBB* JBBBBBBBBB^I 


!!!!!!!!£!K\!!! aaaaaaaa:iaaaaaaaaaaffaBaaaaaaaaa ii a iB 


BBBBBBB BBBBB BB 


bbbbbbbbbbbb 

BBBBBBBBBBBIBBBBflaBBBBB BBaaBBBH^^H 


IBB7BBBBBBBBB tl JBB BB 


BB* BBBBBBB 


IB^1BBBBBBBBBk« BBBBB BBBBB BBBBBBB 


BBB BBBBB BBBBB BBBBB BBBBB BBBBB BBBBB BBJBBB 
BB B BBBBB BBBBB BBBBB BBBBBf Bi B| 


IB«BBBBB 


IBH BBBBB BBBBB B#BBBBB 


BBBBBBBBBB 


■ ■■■BBB' : B BBBBB BBBBB iBBBB BBBBB BlBBBBi 


IBBbB BBBBB BBBBB ’JBBBBI 


llMl^II^III^ 22 ^a^^ 2 a 2 i 222 BBBBBBBBBBBBBBBBBBBBBa ^ BBBBBBBBBB 


20,000 


40,000 
TIME (SEC) 


60,000 


JBBBBI 


80,000 


zr AS A FUNCTION OF TIME 
oe 


-4000 


(BIIIIIBIIIIIBIIIIIIIIIIIIIIIIIIIIIBIIIIIflIIIIIIIIIIBIIIIIBIIIIIIIIIIBIIIIIlj 

Elllllir'MIllllllP'IIIIIIIIBPMIIIIIIIIIlllllllllllllllBIIIIIIIIIIIIIIIII 

Eiiniir iii iiBiiiiiMiiiBiiiiirji iiiimi hiibiiiibf viiiimi r i;miimi 

IBIIIII III ’MIIIMIIllfllllll ilBIIBIBIIIUIIi IBIIIIB JlkHIIIIIMII IIIIIIB 

^Biiiiriiiiiiiiiiiiiift^iiiiirjiiiiiiiiiiJiiL^iiiiiriiiiiiiiiiiiiLiiiiiirj 

hiiiiiiiiiBiiiriiiiituiiiiiiiiiiiiiriiiiitiiiiiJiiiiiiiiiriiiiiLiiiiiH 




-2000 


SllgIBIIBIBIIBilB 


-4000 


iBIHIIBIIMBS 


15*15. 


IIBIIIIBIIHIB 


-6000 


|BMnBimm<mmnSuimm»nnmmmmBmmmmmm^ml 


20,000 


40,000 
TIME (SEC) 


60,000 


80,000 


AS A FUNCTION OF TIME 


FIGURE 16 - 1 UNftR ORBITFR V ORBIT 


- 66 - 



TIME (SEC) 


AS A FUNCTION OF TIME 



TIME (SEC) 


AS A FUNCTION OF TIME 

06J) 


(FT/RAD-SEC) 


-67 - 


•Q. 5 
ro ro 


5000 

4000 

3000 

2000 

1000 

0 

-1000 

-2000 



0 20,000 40,000 60,000 80,000 

TIME (SEC) 


AS A FUNCTION OF TIME 
9M 
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Linear combinations among the Keplerian parameters can be viewed 
as a measure of observability. Hence the slight biases obtained 
in the converged solution can be ascribed to these equal sensi- 
tivites among some state variables and to the aliasing effects 
of the unmodeled short-period variations. 

The second data simulation was generated for a Lunar 
Orbiter III (Apollo type) orbit. This data span contains 
tracking data from the same three stations and is approximately 
10 hours in duration. The epoch date and initial conditions 
for this orbit are: 

Epoch Date : August 30, 1967 20 hours 55 min. 

Initial Conditions : (Selenographic) 

a = 6,457,093 ft. e = .04348376 

I = 20?899211 Q. = 63?970000 

u = 35 4? 05800 M = 194?90793 

The Doppler residuals (Ap) associated with this solution are 
shown in Figure 19. As in the previous solution, these residuals 
are systematic and have the general shape and form of the 
uhmodeled short-period variations. The residuals have a mean 
of -.016 fps and a standard deviation a= .08 f ps . 

The thirteen Keplerian parameters found for this 
solution are the following: 
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e = .04363890 
o 

e ± = -.22143541 x lO -9 

I = 20? 902504 
o 

I I = ,75197331.x 10 -6 deg/sec. 

ft = 6 3? 984837 
o 

= -.14299451 X 10 4 deg/sec 

-10 2 

fi 2 = .52154323 x 10 deg/sec 
to = 354?05553 

to = .23548454 >i 10 4 deg/sec 

-10 2 

to 2 = -. 81426635 * 10 deg/sec 

M = 194?89384 

o 

= .045951285 deg/sec. 

- 1 0 2 

M 2 = .26356956 x 10 xu deg/sec 

The associated semi-major axes found are: 

a Q = 6,457,426 ft. 
a = 6,456,214 ft. 

Variations in the eccentricity and the Euler angles of this orbit 
are shown in Figures 20 and 21. Semi-major axis, mean anomaly, 
and position differences are shown in Figures 22 and 23. As can 
been seen from Figure 20 the inclination parameter has a bias 
error at epoch and a slight slope error. The ascending node 
parameter (see Fig. 21) has only a slight slope error. The semi- 
major axis difference has a bias of 150 ft. The position 
difference has a slight error trend and a bias of 150 ft. 
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MINUS EPOCH VALUE 
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FIGURE 21 - LUNAR ORBITER III ORBIT 



AM (DEG) Aa (FT) 


— 73 — 


•■■■■■■■!■■• .uiaaiki()iaia»\iMMiFiaiMMt ■■■■•■■■raaiii ijiiiiimih 

HIHIPiaHkMllllllllllinHHIir»llllimiillHHIUIIII:INIIIIIIIII 

aiiiiiiiiiiiiiiiaiHii»iitiiisiMiiDiiiiiiiHuiiiaiiMiiiMniai»iiiii 

nKiifiiifiiiiniiiiiiiiiiuiaiiHiiimibiarniMiikiiHiiuiiiiiiiii 

aiiiiiiiaiiaiiiim»iiiaiFiiiiii)iaiijiiiiimiiMiaiiaiHia?iaiiaiaM 

iiaaiiiiaifiii'iiiiiiiKiiiiiaaaiiii'iaaaaii'iiiniiniiiiiiiiiiiiiiiiiaii 

Maaiaiiaiiaiaiiaiaaaiaaai»ainiaaaaiaaiauaiaaia?iaauaa(iaaKaaaiiaaii 

aniiniiiMiiEiaiiaiiaiaiMiiiiiiNiiLiinianiiiiiaiiaiiiiiiiainiiia 

miaaaiiMniniami)iHiaa?aaanaaaaiaaiiimuianiiiraiiMiiaaaiia 

aiMiiKaiaii»iaiFiiiiiia«iHilaiiiiiiiiHiii?iiiyiiiiiiiiiiiiiiiiiiiii 

aiMiiaiiiJiaaaiaiiaiiaiaiiaairiia^iaaiaakaiaaaiaaBaaaaaaiiaaaaaaKaaaav 

BBBBBIBBBkBBBBBBBBIBBBBBBBBBBBBBABIBB.BBB^BUBBBBBBBfeBBBBBJBBBBBBRBBnBBBBB 
1BriBBBBBB^BBBBBbBBVBBBBI.'B^BBIBBKBbBBBBBBBaBBBBB'JB?BBBBBBIiflBBBflBBBIBaBBr 
BVlBBBBBBtlBIBBBriBBBBBBBflBBBBBSBffBBBBflFkBBBBBBiPBBaBBMbrBBflBataflBIBM 
artBBBBBi IBBiBBf ILBBBtBl' lBBBBBLB^BBlBBf'sBBBBBBKBJBBBBBt 7BBBBBBia*BBBBBB 
^BaaaBaBaaBBBaB^^BBBBBBiBBBBBB.KfllBBBBBBrfBBBBaBbC'lBBBBBbiBBaaagfbaBBBBB 
■LBBBBBBBBBBaaBBkiBBBBBBaBBBBBBBLBBflBBBBBBBBflBBB JBBBBBBBBBBBBBt jBBBBBB 
[■■■■■■■■■■■BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB 


24,000 


TIME (SEC) 


DIFFERENCE IN SEMI-MAJOR AXIS AS A FUNCTION OF TIME 


I l I I l i l I i i i l i i i J i i i l M l l M i l l i l I ll l I I 


lllEMlliiBiMS!MS2!iSi!!i!S!!5! B M£ BBBBBBBBBBBBBLBBBBBBBBBBBBBBB ' ? i BBBBBI 

nS:L!SS!SS!!!! BB !!! BBBBBB !!! B !! Bfc !f; BBBBBBBBBBBBB - BB « BBB « BBBBBBBB -« BB * 

SSiwMliBgilBaBjaaa^lBpBpBB'lj^MJBJ^iaaBBBBBBBBBBBBSBBBBBBBBBBBBBBBJPflBBBM 


iJBLBaaBBaBBBBaaLBBaBBBBaBBBBaaB^BB:iBaBBaBBBBBBa*BLaaBBBaaBBBBBrBB*BBB 
■MaaaaaBaaiaaaflB^BBiaaaaaaaBaaaBiaflaaBaaaBBaBaflEjaaBBaBaaBBBBBBraa?aBa 
IBTBBL1BBBBBBBBBBB1B BBBB BBBB BBBBBBBBB^BBBBBBBBBB BpBBl JBB BBBB BBBBBNBBfcBBB 
8fi*BBBBBBBBBBBBB.BBBBBF > BBBBflBBBBB’BthBBBKBBBBBBBBBBaTBBBBBBBBBBBBBBB4BBBBBB 
•■^BBBlBBBBBBBBBB^BBBBBBBBBaBBBBBtvBBBBBBBBBBBBBB B.IBBBnBBBBBBBBBBB JBBL1BB 
■rBBBBBaaBBBaflBBa^. BBaB^aBBBBBBBBB dBBflLflBaBBBBBBBI JBB BBBB BBBB BBBB^ BBBBBBB 
BrBBBB?BBflBBB8BBriBBBBBBBBBBBBBB r BBBBBBBBBBBBBBB BBBBFBBBBBBBBBB^/BBB BIBB 
f iBBBBBBBBBBBflBB'- BBBBBBBBBBB BBBF JBBBBI JBBBBBBBB' iBBBBUBBBBBBBBB' . BBBBBBBB 
BBBBBPflBBBBBBBf JBBBBBy B BBBB BBB BBBBBBBBBBBBBBii BBBBBBBBBBBBBBrj BBBBB7BB 
JBBB BBBB BBBBBBB VBBBBBBPBBBBBBBB 4 BBBBBB ‘IB BBBB BBF BBBBB BBBBBBB BBB>lBBBBBBBBB 
•BBB BBV ^ BBan IlkiBBB I ■ BBT B BBBBBBBB BBBBBBBB I BBBIBBBBBIBBBJBBBBBBrBB 
IBBB BBBB BBBB BBPB BBBB BBr VBBBB BBBBBBBB BBwBBBBB BBBBBBBB BBBBBBBB Bt'JB BBBBBBBB 
BBBBBBB JBBBB BBBBBBBB BBBBBBBB BBBBBBBBBBrBBBBBBBJB BBBBBBBB BBBBBrBBBBUBBBK'JB 
■BBBBBBrBBBBBBLBBBBBBBB^BBBB BB ’JBBBBB BBBBBBBB BriBBBBB BBTB BBBB BBBB BBBBBBBB 
BBBBBBBr BBBBBBBB BBBBBBBr/BBBB BBBBBBBB BBBBBBBB BBBBBBBB BBM BBBBBBBB BBBBBBBB 
■BBBBBBB -^BBBBBBBBBBflBBBk BBBBBBBBBBBBBBB i BBBBBLlBBBBBBBBB^BBBBBBBBBBBBBBl'B 
■BBflflflBBk IBB BUB BBBBBBBB i BBBBrBBBBBBBBBB ’/BBBBBBB BBBBBBBk 1BBB BBBBBBBB BBB, . 
■BBBBBBBBOBBBBBBBaBBBBBBB 1BBBBBB BBBBBBBB L*BBBBJBB BBBBBBBB ^BBBBBBB BBBBBBBB 4 
■BBBBBBB BcBB BkJB BBBBBBBB B^BB a:* BB BBBBBBBB LJBBBBBflaBBB BBBI ^BBBBBBBBBBBBBBB^ 

■ BBB BBBB BB 'IB BtaBB BBBB ■ BBB B laBB BUB BaBBBBB.BB BV BBrVBB BBBBBBBB BJBBBBBB BBBBBBBB 
BBBB BBBBBBBB BBBBBBBB ■BBBfllbBBBBBBB,flB.B.B.BBBBBBBBBBBBBBBBBBBB*]BB>^BBBBBBBBBBB 
BBBBBBBB BBBBBBBB BBBBBBBB BB -IB BBBBBBBB BBBBBBBB^BBBaBBBBBBBBPBBBBBBBBBBBBBB 
IBBBBBBBBBcB t IBB BBBBBBBB BB^B TBBB BBBBBBB BBBJBBBBBaBBBBBBBBBBBBBBBBBBBBaBB 
■BBBBBBBBBrBPBBBBBBBBBBBBBirtBnBflBBBBJBBBBBBB>lB*rfBBBBBBBBBBBBBflrBBBBBBBBBBBB 
■flBBBBBBBBBBNBaBBBaBBBBBBBk’BHBBaBBBBBaaaBBy'B JBBB BBBBBBBB BBJrBBBB BBBBBBBB 

■ BflflBBBflBBB:! ^IBB BBBB BBBB BBfc* UBBBBBBBBBBB BBC* V BBBB BBBBBBBB BB/VBBBB BBBBBBBB 
iBBBBBBB BBBBBBB^ BBBBB.B.BBBBBBBBB- iBBBBBBBB BBBBBB kJBBBB BBBBBBBB 


12,000 


24,000 


36,000 


TIME (SEC) 


DIFFERENCE IN MEAN ANOMALY AS A FUNCTION OF TIME 


FIGURE 22 - LUNAR ORBITER III ORBIT 
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DIFFERENCE IN POSITION AS A FUNCTION OF TIME 


FIGURE 23 - LUNAR ORBITER III ORBIT 


75 


Again an examination of the functional properties 
of the partial derivatives used in this solution reveals that 
three of the state variables have virtually identical sensi- 
tivities to the data (see Fig. 24-26) 


3 P = 3 p = 3 p 

3^ Tto 3M 


(3.52) 


. T -1 

High correlations in the solution [H w H] matrix were found 

to exist between the following sets of parameters: 


Correlation Coefficient Parameter Pair 
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Again the slight biases and trend errors experienced 
in the convergence are a reflection of the equal sensitivities 
among the state parameters and the aliasing effects of short 
period variations. 
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FIGURE 24 - LUNAR ORBITER III ORBIT 
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FIGURE 26 - LUNAR ORBITER III ORBIT 


79 


The overall results of these simulations show that 
the empirical orbit determination method can be used to 
accurately estimate long-period orbital element variations 
arising from the non-central features of the moon, from earth 
and sun perturbations, and from solar radiation. This analysis 
shows the largest single error source is the lack of separability 
of dynamical effects from the Doppler measurements. Analysis 
has shown (for Doppler data) the only way of enhancing these 
sensitivities is to choose earth based tracking stations with 
the greatest north-south separation relative to the lunar orbit 
and to have simultaneous tracking coverage from these stations 
whenever possible. 

The set of Kepler element rates, orbital elements, and 
associated weighting matrices obtained from this pseudo data 
analysis will be used to determine a set of triaxial lunar 
coefficients. The results of that procedure will be discussed 
in the next chapter. 



CHAPTER IV 


HARMONIC ESTIMATION 


The gravity field determination is performed in 
a second weighted least-squares processor which uses as input 
the Kepler element rates, the estimated Kepler elements, and 
a weighting matrix and outputs a set of spherical harmonic 
coefficients (see Figure 27) . The perturbation equations are 
of the form 


k = f (k,p,t) 


(4.1) 


where p is the vector (nxl) of gravity coefficients. 



(4.2) 


The error function to be minimized for this estimator is: 


e = (k - k*) T A(k - k*) (4.3) 

where A is the (5x5) weighting matrix introduced in the pre- 
vious chapter, k* is an estimate of k. 

LEAST SQUARES PROCESSOR 

Since the gravity parameters appear as linear 
functions in the perturbation equations, then (4.1) can be 
expressed in the following form: 
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FIGURE 27 - GRAVITY COEFFICIENT ESTIMATION 
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k* = F(k)p* 


(4.4) 


where F is a (5xn) matrix of partial derivatives of the Kepler 
element rates with respect to the gravity coefficients: 


F (k ) = 



3e 3e . . . 3e 
9C £m ' 9S 21 9S lm 


(4.5) 


3M 

3C 


20 


3M 

3C 


.fin 


3M 
' 3S 


21 


3M 

3S 


-fin 


The error function can be written in terms of the gravity 
coefficients as follows: 


e = (k - Fp*) T A(k - Fp*) (4.6) 

The minimum of the error function is found by setting 6e = 0. 
The resulting equations are: 



T — 1 
[F AF] 


(4.7) 


(nxl) (nxn) 


(nxl) 
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where p is the best estimate of the lunar gravity parameters. 
TRIAXIAL FIELD DETERMINATION (PSEUDO DATA) 

In order to demonstrate the dynamical continuity 
between the empirical orbit determination and harmonic estimating 
processors, the two pseudo data convergences presented in the 
previous chapter will be used to find the triaxial (C 2 g and C 22 ) 
lunar coefficients assumed in the source trajectory. The 
nominal values of these terms used for pseudo data generation 
are as follows: 


C 2Q = -2.07108x10 
C 22 = .20716x10 


> (4.8) 

) 


Since the orbital period of the Orbiter V satellite 
is 3.2 hours and the data span 21 hours, this converged solution 
contributes seven sets of Kepler element rates and elements 
to the harmonic estimator. The Orbiter III satellite has an 
orbital period of 2.1 hours and a data span of 10 hours, con- 
sequently this solution contributes five sets of Kepler elements 
and element rates. 

The numerical values for the lunar triaxial coefficients 
as obtained from the pseudo data Orbiter III and Orbiter V con- 
vergences when used as input to the harmonic estimating processor 
are as follows: 


C 

C 


20 

22 


-2 


.09x10 


-4 

-4 


.209x10 


(4.9) 
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Since this entire process is one of parameter identifica- 
tion (See Appendix E) and many of the state variables are subject 

to systematic errors and biases, it is questionable whether 

T -1 

the terms on the diagonal of the [F AF] matrix can be 
regarded as variances of the C 2 Q and C 22 estimates. The 
principal detriment to making this assumption is the fact 
that A weighting matrix used in the process is not truly the 

covariance matrix of the input Kepler element rates. The 

. . T -1 

normalized non-diagonal terms m the [F AF] matrix do reflect 

the true measure of correlation between the C 2 Q and C 22 terms. 

The correlation coefficient between these two terms for this 
harmonic determination is p = .08. 

Analysis shows when each of these solutions was used 
separately for a coefficient determination the numerical re- 
sults obtained varied. The most significant aspect of using 
only single satellite solutions for gravity determination is 
that the correlation coefficient between the parameters in- 
creases. The harmonic estimation results from using each 
satellite independently are given below: 


SATELLITE 

ESTIMATED 

CORRELATION 

NUMBER 

HARMONICS 

COEFFICIENT 

ORBITER V 

C 20 = -2.056xl0~ 4 
C 22 = .219 x10 _4 

p= . 98 

orbiter III 

C 20 = -2.044xl0“ 4 
C 22 = .668xl0 -5 

p = . 99 
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The small correlation coefficient associated with 
the multi-satellite solution has an important physical inter- 
pretation. Basically a high correlation reflects the in- 
separability of gravity effects in the harmonic estimator. 

Since the gravity coefficients each affect different orbits 
in different ways, using solutions from as many orbits as 
possible reduces this inseparability. Optimum separation of 
dynamical effects is achieved by using data from orbits of 
many different inclinations. The gravity estimation errors 
associated with multi-orbit solutions are then largely a function 
of the degree of accuracy of the estimated Keplerian states 

and rates. For the case of single orbit solutions, the 

T -1 

presence of high correlations in the [F AF] matrix also tends 
to confuse the estimation process. This notion of separability 
of dynamical effects becomes very important in the actual deter- 
mination of the lunar field. 

In order to obtain a quantitative measure of the tri- 
axial solution given by (4.9) a comparison is made between it 
and the nominal triaxial values (4.8) using numerically inte- 
grated eauations of motion. The comparisons, covering a 
one day period, are made both for the Lunar Orbiter III and V 
orbits. Differences in each of the Kepler elements, position, 
and velocity are shown in Figures (28-35) . The position errors 
developed over the Lunar Orbiter V trajectory are small and 
attain a maximum value of 60 ft. The velocity errors are also 
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FIGURE 28 - LUNAR ORBITER V ORBIT 
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INCLINATION DIFFERENCE VS. TIME 



ASCENDING NODE DIFFERENCE! VS. TIME 


FIGURE 29 - LUNAR ORBITER V ORBIT 
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FIGURE 30 - LUNAR ORBITER V ORBIT 
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FIGURE 32 - LUNAR ORBITER III ORBIT 
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FIGURE 33 - LUNAR ORBITER III ORBIT 
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FIGURE 34 - LUNAR ORBITER III ORBIT 
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small, reaching a peak value of .05 fps. For the case of the 
Lunar Orbiter III trajectory the position errors are even 
smaller, attaining a peak value of 25 ft. and the velocity 
errors attaining a value of .025 fps. 

DISCUSSION 

The dynamical compatibility of this empirical 
selenodesy method has been established in this chapter. It 
is significant to note at this point that the method is de- 
pendent on the following: 

1. Doppler tracking data of free-f light quality 
(thrust-free) . 

2. Tracking coverage relative to the satellite 
orbit which provides the best geometrical con- 
figuration possible (good north-south station 
separation) . 

3. Tracking of satellites from many different 
inclinations so as to attain the best overall 
coverage of the moon. 

These points are re-emphasized since in the next chapter of 
this investigation the empirical method is applied to actual 
Doppler data. Of the three points mentioned, none are actually 


achieved for the case of real data. 



CHAPTER V 


DATA ANALYSIS 


In the previous two chapters, the theory and equations 
for this empirical selenodesy method have been developed 
and a controlled pseudo data analysis has been presented to 
illustrate the dynamical consistency of the method. The purpose 
of this chapter is to present and discuss the results obtained 
when this method was applied to actual Doppler tracking data. 

DATA SET UTILIZED 

As mentioned earlier in this study, large amounts of 
Doppler tracking data were acquired during the lunar orbits of 
both the Lunar Orbiters (I-V) and the Apollo (8, 10, 11, 12, 14, 
and 15) missions. Almost all the tracking data acquired during 
the photographic portions of the Lunar Orbiter missions includes 
propulsive attitude control maneuvers performed at such a high 
frequency (about every three hours) that these data cannot be 
used for selenodesy purposes. Even the Lunar Orbiter data from 
the extended mission phases (primarily that used in this analysis) 
contains some minor propulsive thrusts. Data from all the 
Apollo 8 mission and large portions of the Apollo 10, 11, 
and 15 missions contain propulsive thrusting. The extended 
mission phase Lunar Orbiter data was used to determine the 
lunar gravity field presented in this study since it is not 
only the largest but also the most complete (but far from 
complete in an absolute sense) data set gathered to date. 
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The free-f light quality data from parts of the Apollo 11, 

12 and 14 missions were used as control data to test the 
quality of the lunar field obtained. The small amount of 
Apollo data which is of free-flight quality cannot be used in 
the method since these data are from either nearly circular, or 
nearly equatorial orbits and the method as it has been developed 
becomes near singular for orbits of this type. 

The epoch times, length, and number of tracking stations 
of the various data arcs used in this analysis are listed in 
Table III. No data from the Orbiter IV satellite were included; 
since earth perturbations were the most dominant for this orbit, 
it was assumed to have only minimal lunar gravity informa- 
tion content. The length of the data arcs used in orbit deter- 
mination solutions varies from a minimum of about eight hours 
to a maximum of thirty -six hours. 

ANALYSIS OF ORBIT DETERMINATION SOLUTIONS 

Two different sized parameter sets representing the 
satellite Keplerian state in lunar orbit were used. The first 
set contains eleven parameters (^ e Q r, e^, I Q , 1^ 

M q , an< ^ was use< ^ exclusively when the Doppler data span was 

twelve hours or less. The second solution set contains thirteen 

parameters / ^1 ^ ^2 r ^o * ^"1 ^ ^o f t ^2 * ^o 9 ^1 ^ ^2 * ^o r ^1 ^ ^2^ 
and was used when the data span was greater than twelve hours. 

Analysis with pseudo data showed that this choice of parameters 
should be adequate to model the long-period variations of the 
satellite . 
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SATELLITE 

EPOCH (DAY, 

MO. , 1 

m. ) 

LENGTH (SEC) 

ORBITER I 

31.41180556 

Aug . 

1966 

75,989 


1.32430556 

Sept. 

1966 

88,889 


4.22916667 

Sept . 

1966 

84,000 


13.84722222 

Sept . 

1966 

81,000 


14.81597222 

Sept . 

1966 

81,000 


15.78263889 

Sept. 

1966 

108,670 

ORBITER II 

8.85938421 

Dec. 

1966 

92,700 


10.55208333 

Dec . 

1966 

32,760 


29.61805556 

June 

1 Q 67 

34,260 

ORBITER III 

20.60000000 

Feb. 

1967 

93,290 


21.68055556 

Feb. 

1967 

86,000 


27.97222222 

Feb . 

1967 

89,280 


2.38194444 

Mar . 

1967 

75,000 


5.16666667 

Mar . 

1967 

110,800 


6.83958333 

Mar. 

1967 

97,740 


24.21527778 

Mar . 

1967 

31,380 


11.48611111 

April 

1967 

46,080 

ORBITER V 

18.31111111 

Aug . 

1967 

63,000 ! 


19.09027778 

Aug . 

1967 

88,000 


20.11111111 

Aug. 

1967 

119,450 


21.49652778 

Aug. 

1967 

119,068 


24.21180556 

Aug . 

1967 

87,928 


25.22847222 

Aug . 

1967 

59,340 


26.25208333 

Aug . 

1967 

85,000 


27.28819444 

Aug. 

1967 

70,468 


2.52777778 

Oct. 

1967 

29,880 


3. 54166667 

Oct. 

1967 

43,200 


17.38194444 

Nov. 

1967 

38,040 


21.31944444 

Nov. 

1967 

33,568 


29.86111111 

Jan . 

1968 

112,618 




roromrororol co ro cm! rorororororocNCM I rorororororororocNCNCNCNro 
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TABLE IV PHYSICAL CONSTANTS 


I. ASTRODYNAMIC CONSTANTS 


Gravitational Parameters: 

p^ = .1731300417087798 x 10 15 Ft 3 /sec 2 
p Q = .1407646853278542 x 10 17 Ft 3 /sec 2 
p = .4686697671960888 x 10 22 Ft 3 /sec 2 


Mean Lunar Radius : 

R = .570239501312336 x 10 7 Ft. 
e 

Angular Velocity of Moon's Rotation: 

y = .2661703316891657 x 10" 5 Rad/sec 
Velocity of Light in a Vacuum: 

c = .9835711942257218 x 10 9 Ft/sec 


II. STATION LOCATIONS 


STATION | 

GEOCENTRIC COORDINATES 


Radius, Ft 

. Latitude, Deg. 

Longitude , Deg. 

Goldstone, Calf. 
(DSS12) 

20,905,479 

35.118640 

243.19483 

Woomera , Australia 
(DSS41) 

20,907,326 

-31.211390 

136 .88779 

Madrid, Spain 
(DSS61) 

20,898,911 

40 .238540 

355.75129 

Madrid, Spain 
(DSS62) 

20,898,927 

40 . 263490 

355.63246 


III. DOPPLER DATA WEIGHT W = -t , a = .00213 Ft/sec 

a 

IV. SOLAR ACCELERATION F = 3.973 x IQ -7 Ft/sec 2 
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Table IV contains a list of the tracking station 

locations^ and all the astrodynamical constants^ assumed known 

and fixed in the orbit determination program. The positions 

and relative velocities of the earth, sun, and moon are obtained 

32 

from Jet Propulsion Laboratory Ephemeris Tape DE-19. The 
least-squares orbit determination program reads the necessary 
information from this tape during the convergence process. 

The rate of convergence for each orbit determination 
solution varied as a function of the orbit, the data arc length, 
and the number and locations of earth based stations used. In 
general each convergence process took about five to six iterations. 
The numerical particulars associated with a typical solution for 
each of the Lunar Orbiters (I, II, III, and V) will now be 
discussed. 

Representative residuals from Lunar Orbiter I, 

II, III, and V orbit determination solutions are given in Figures 
36 and 37. The residuals associated with each of these convergences 
have the statistical properties listed in Table V. 

TABLE V 


SATELLITE 

RESIDUAL MEAN (fps) 

RESIDUAL STD. DEVIATION (fps) 

ORBITER I 

-.010 

.108 

ORBITER II 

-.012 

.090 

ORBITER III 

.006 

.053 

ORBITER V 

.003 

.068 


As can be seen from the figures the residuals associated with 
the Lunar Orbiter I and II convergences presented have a larger 
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peak amplitude. Each of these residual plots also possesses 
points of irregularly large amplitude relative to the remainder of 
the span. These irregularities in the residuals most likely 
correspond to low thrust attitude control maneuvers performed 
by the spacecraft. 

. T -1 

An analysis of the [H WH] matrix for each of these 
solutions shows that many of the Keplerian parameters are 
extremely highly correlated. Table VI presents a summary of the 
highest correlations for each of these four convergences. 


TABLE VI SOLUTION CORRELATIONS 


| SATELLITE 

| PARAMETER PAIR 

CORRELATION COEF . 

i 

ORBITER I 

j 

Q o' 

00 

o 

-.9999 | 


r 

“1 

.9999 | 

ORBITER II 

V 

n ! 

o 

-.9908 ] 


a o' 

00 

O 

-.9998 

i 

9 

u i 

-.9995 

ORBITER III j 

V 

Q 

o 

-.9938 

1 

1 

V 

00 

o 

-.9946 > 


n o ' 

00 

o 

-.9998 | 


: fl l r 

W 1 

-.9985 ; 

ORBITER V 

V 

0 

H 

a 

o 

.9948 


V 

00 

o 

-.9903 j 

j 

SIq, 

00 

o 

-.9843 : 


9 

h' 

.9948 ! 

s 


These correlations represent the largest found in 
T -1 

each [H WH] matrix. There were other correlations present of 
large magnitude (p = | . 9 J ) . For the case of each Lunar Orbiter 
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satellite presented almost the identical correlation pairs 
reoccur. These correlations reflect the basic difficulty of 
using Doppler measurements at lunar distances to separate the 
dynamical properties of the Euler angles of the orbit. It should 
be noted that this situation is particularly amplified for the 
case of the Lunar Orbiters since in most cases the satellite was 
only tracked by one earth-based station, thus losing the geome- 
trical enhancement of a second or third tracker. 

The thirty data arcs used for orbit determination 
solutions contributed one hundred ninety-nine sets of states and 
rates for lunar gravity determination. The factors used in defin- 
ing the degree and order of the lunar field to be determined from 
these sets of Keplerian rates and elements will now be discussed. 

A PRIORI COEFFICIENT SELECTION 

33 

The work of Muller and Sjogren has provxded conclusive 
evidence that rather large near-surface mass concentrations 
("mascons'j are present in the near-side lunar maria regions. 

The existence of these non-central mass concentrations has a sig- 
nificant impact on the application of Equation (2.20) for 
describing lunar gravity. Accurately representing a "mascon" 
moon would require spherical harmonic coefficients of 
high order and degree. In mathematical terminology, the presence 
of "mascons" causes the convergence rate of (2.20) to be very 
slow. 

In theory then the proper approach to modeling the 
lunar gravity field is to seek a spherical harmonic coefficient 
set large enough in degree and order to represent all the non- 
central lunar features. In practice however, due to the 
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incomplete Doppler data set and due to a lack of far side 

Doppler measurements, lunar gravity solutions involving large 

numbers of harmonic coefficients have high correlations in 
T -l 

the [F AF] matrix and in general have poor overall numerical 
characteristics . 

Analyses were made using the harmonic estimating program 
with pseudo data input from numerically integrated long-period 
perturbation equations (assuming a nominal seventh degree and 
order field). Long-period trajectories for each of the Lunar 
Orbiter missions were simulated. When a solution set of degree 
and order seven or larger was sought, the nominal values of 

the field (as could be expected) were recovered with good accuracy. 

. . T -1 

An analysis of the correlations in the [F AF] matrices for 

these solutions (for Lunar Orbiter tracking coverage) showed 

that a great many harmonic coefficient pairs (C 20 & ^40' 

C 30 & ^50 etc ) were very highly correlated. This correlation 

is totally a consequence of the incomplete data set used. 

When subset gravity solutions (for example a complete fifth 

degree and order field) were sought from the Keplerian rate data 

generated from the seventh degree and order field, the numerical 

values obtained were very different from their nominal values. 

Basically, the higher degree harmonics which had been omitted 

from the solution set were aliased into the lower ones due to the 

existing large correlations. Had a complete data set (data 

covering all latitudes and longitudes) been used, then ortho- 

. T -1 

gonaiity would be induced m the [F AF] matrix and the 
subset values recovered would be the nominal ones. Since the 
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spherical harmonic expansion (2.20) is essentially a three 

dimensional Fourier series in spherical coordinates, this 

34 

orthogonality can only be achieved when: 


if or m^k 


Y hk d ° 


(5.1) 


4 tt ( Z+m ) ! 


sphere 


^(l-m)I (21+1) (2-6 Qm ) 


for t=h and m=k 


where: Y^ m (<p,X) = P^ m (sin<j) ) [C^ m cos mX + S^ m sin mX ] 

and where the Kronecker delta 6 is equal to 1 for m=0 and 0 

om ^ 

for m^O. 

This discussion of the "mascons" and the lack of 

orthogonality associated with the existing Doppler data is 

introduced as background for the rational process used in choosing 

a harmonic coefficient solution set. The basic strategy 

assumed in this analysis is to obtain the largest coefficient 

set possible while incurring a minimum of high correlations in 
T — I 

the [F AF] matrix. It is a foregone conclusion that, with 
the data available at this time, it is not possible to determine 
a lunar gravity field which truly represents all the localized 
fine structure near-surface mass inhomogeneities. The only 
attainable goal of this data analysis then is to derive a more 
accurate global lunar gravity model. 
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A LUNAR GRAVITY FIELD 

Solutions varying from degree and order three to 
degree and order seven were attempted from the one hundred 
ninty-nine sets of Kepler element rates. All solutions ob- 
tained above degree and order four contained very large 

T -1 

numbers of parameter correlations in the solution [F AF] 
matrices. When these solutions were applied to the tracking 
data from the Apollo orbits, both the fit and prediction 
characteristics obtained were very poor. As a result of 
this situation, the lunar harmonic coefficient set determined 
in this study is of order and degree four. 

Analysis of the numerical characteristics of the 
full fourth degree and order solution revealed two important 
points. First, the C^q and C^q zonal coefficients were still 
highly correlated (p = .86). Second, the solution contained 
very little direct information in determining the C 22 harmon- 
ic. In performing the least-squares gravity estimation, the 
T — 

entry in the F Ak vector associated with C 22 was essentially 
zero (other components were of significantly larger magni- 
tudes) . Hence the estimate of C 22 was dominated by correla- 

T - 1 

tions present in the [F AF] matrix. 

In order to circumvent these numerical problems , 

both the C 2 q and C 22 terms were fixed in the gravity deter- 

. 35 

mination to values obtained by Koziel in studying the physi- 
cal librations of the moon. The values used are the following: 
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The fourth order and degree field lunar gravity 
field (with C 2 Q and fixed) is given in Table VII. The cor- 

relation matrix associated with this solution is given in Table 
VIII. As in the pseudo data analysis previously presented, 

this gravity estimation process is subject to systematic 

T -1 

errors and biases. Hence the diagonal terms of the [F AF] 
matrix cannot be regarded as variances of the estimated terms. 

The residuals associated with each of the five element 
rates and this fourth order field are shown in Figures 38-40. 
These residuals have the statistical properties listed in 
Table IX. 


TABLE VII 

GRAVITATIONAL FIELD DETERMINED FROM LUNAR 
ORBITERS I, II, III, AND V 


£ 

m 

4 

C „ x 10 * 
£m 

4 

S„ x 10 
£m 

2 

0 

-2.0560* 

— 


1 

.0537 

.0617 


2 

.2258* 

-.0017 

3 

0 

.2216 

— 


1 

.3575 

.0820 


2 

.0210 

.0340 


3 

.0301 

.0055 

4 

0 

.0543 

— 


1 

-.0677 

.1195 


2 

.0443 

.0106 j 


3 

.0136 

.0066 ! 

1 

4 

.0027 

.0043 j 


*Fixed in solution. 
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TABLE IX 


Kepler Element 
Rate 

Residual Mean 

Residual Std. 
Deviation 


-9 

-8 f 

e 

-.272 x 10 

.345 x 10 

• 

i 

-.836 x 10 7 Deg/Sec 

.170 x 10 ^ Deg/Sec 

n 

.127 x 10 5 Deg/Sec 

.265 x io 5 Deg/Sec j 

• 

0) 

-.965 x io 6 Deg/Sec 

.374 x io Deg/Sec \ 

\ 

M 

— 7 

! -.915 x io Deg/Sec 

i 

— 1 

.327 x 10 Deg/Sec j 


As can be seen from the Kepler element rate residual plots, the 
errors are systematic in each case. 

Equipotential surfaces have been calculated for this 
lunar gravity field and are shown in Figures 41 and 42. These 
surfaces are computed by finding the radial deviations to a spher- 
ical potential (generated with the field point at the mean lunar 
radius) . The variations are quantized in thousand foot incre- 
ments. The basic equipotential surfaces of this gravity field 
are those of a triaxial ellipsoid. The solid line on the 
surfaces indicates the equipotential line for the reference 
potential (zero deviation from spherical potential) . These 
surfaces show three large areas of potential excess. The first 
of these is centered about latitude <|> = 25°N and longitude 
X = 10 °E. This region very closely corresponds to the Mare 

Serenitatis region of the moon. The two other areas of poten- 
tial excess are located at latitude (j> = 5 °S, longitude X = 117 °E 
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and latitude <t> = 10°S, longitude X = 170°W respectively. Neither 
of these two areas corresponds to an identified lunar maria 
region. 

The second-degree harmonics in the potential are direct- 
ly related to the moments and products of inertia of the moon. 

The relations between the gravity coefficients and the moments 

3 6 

and products of inertia are: 


C 


20 


m d R e 


A+B 

2 



C 


21 






22 ' n * e 2 


(B-A) 


22 2 

m d R e 


J 


(5.3) 


where A, B, C are the three principal moments of inertia 

(A = If B = I , C = I ) and D, E, F are the products of 
xx yy z z 

inertia (D = I . E = I , F = I ) . Since equation 5.3 con- 
y ^ x z xy 

tains five equations in six unknowns (A, B ,C , D ,E , F) one addi- 
tional relationship is needed. From studies of the lunar 
physical librations the quantity 

C—A 

6 = ^ (5.4) 

has been determined. The numerical value used is that computed 

37 -4 

by Koziel from heliometer observations (B = 6.294 x 10 ) . 

Hence given the five second degree harmonics and 3, the follow- 
ing set of principal moments and products of inertia are found: 
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D = .617 x 10 5 m.R 2 

^ e 

E = .537 x 10 -5 m. R 2 

V e 

F = -.17 x 10 -6 m.R 2 

^ e 

The imprecision in the harmonics and the simplifica- 
tions in the theory relating 6 to the inertias make the quality 
of these numbers somewhat poor. 

EXTRAPOLATIONS 

In order to measure the orbit determination and 
prediction qualities of the field obtained, it was applied 
to Doppler data not used in its generation. Specifically, 
the data used was acquired during the orbits of the Apollo 11, 
12, and 14 missions (Apollo 15 data is not available at this 
time) . The orbit determinations were performed using a standard 
least-squares processor which obtains a best estimate of a 
rectangular state vector at some epoch. 

X q = [H T WH] -1 H T W Ap (5.5) 

This process is identical to that given by (3.31) with the 
exception that is a six vector of position and velocity 
and the H matrix contains the partial derivatives of the 
Doppler with respect to the rectangular state at epoch. This 
best estimate of the state is then used to predict the Doppler 
outside the span of data used for the convergence. 

Orbit determination solutions were obtained by fitting 
one front side pass of Doppler data from several stations. This 


A = .3983208 m^ 2 

B = .39841118 i\R e 2 

C = .39857195 m,R 2 
K e 
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particular data length was chosen since it puts maximum stress 
on the orbital period prediction capability of the model. Once 
a converged solution is obtained, the Doppler data is predicted 
for the next three orbital periods to test the extrapolation 
capabilities of the model. 

The data used from the Apollo 11 and 12 flights are 
from near-circular orbits with radius vectors of about 60 
nautical miles (n.m.) above the lunar surface. The Apollo 14 
data is from a slightly elliptical orbit (e = .0258) with an 
apolune of 60 n.m. above the surface and a perilune of 8 n.m. 
above the surface. The data from the Apollo missions, since 
they are collected from orbits very near the lunar surface, 
reveal many gravitational perturbations not present in the 
Lunar Orbiter data. 

Converged solutions were obtained using the fourth 
degree and order gravity field. Doppler residuals for both 
the one pass fit and three passes of prediction for each of 
these Apollo orbits are shown in Figure 43. The Doppler 
residuals in each of these convergences are still systematic 
and an order of magnitude above the noise level of the data. 

The three orbital period prediction is characterized by a 
secular growth (period errors) in the residuals for each case. 
Both the systematic nature of the regressed Doppler and the 
growth in Doppler residuals during prediction reflect the in- 
complete nature of this fourth order field. 


In order to obtain a relative perspective on the 
quality of orbit determination and prediction attainable, 



Ap(FPS) Ap(FPS) Ap(FPS) 
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FIGURE 43 - FOURTH DEGREE AND ORDER LUNAR GRAVITY FIELD 
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convergences were also performed on these data using two other 

gravity models. The first model, the Ll field (see Table II)> 

is used by the Manned Spacecraft Center for Apollo mission 

3 8 

planning. The second model, developed by Liu and Laing, is 
a fifteenth degree zonal and eighth order tesseral model (84 
harmonics). This model represents the latest result from the 
indirect (long-period) analysis method. 

The orbit determinations were again performed by 
fitting one pass of data and predicting the Doppler for the 
next three periods. The residuals associated with each of 
the fit and predictions for these two gravity fields are shown 
in Figures 44 and 45. Again both the Doppler residuals for 
the fit and prediction, for each field, have systematic errors. 
Both of these models, especially the fifteenth degree field, 
are characterized by secular growth in Doppler residuals of 
the prediction. Table X lists the statistical properties of the 
one pass Doppler residuals for all three models. 

TABLE X 


ORBIT 

MODEL 

RESIDUAL MEAN (FPS) 

STD. DEVIATION (FPS) 

Apollo 11 

4x4* 

-.011 

.166 


Ll 

-.0015 

.039 


15 x 8 + 

-.023 

.463 

Apollo 12 

4 x 4 

-.0029 

.136 


Ll 

-.0042 

.104 


15 x 8 

-.022 

.353 

Apollo 14 

4 x 4 

-.0016 

.187 


Ll 

-.0055 

.159 


15 x 8 

-.0304 

.526 


*Fourth degree and order model 
tFifteenth degree and eigth order model 
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For these Apollo orbits, the LI field achieves slightly better 
convergence statistics than the fourth degree field. However, 
the fourth degree field predicts the Doppler with an error rate 
of about 50% less than LI. Of the three fields compared, the 
fourth degree field most accurately describes the orbital varia- 
tions arising from the non-central features of the moon. 



CHAPTER VI 


SUMMARY AND CONCLUSIONS 


This study has presented an empirical method for 
determining the spherical harmonic coefficients of the lunar 
gravity field. This method uses a two-step data reduction and 
estimation process. In the first step, a weighted least-squares 
empirical orbit determination process is applied to Doppler 
tracking data to estimate long-period Kepler orbital elements 
and rates. Each of the Kepler elements is represented by an 
independent time function. The long-period perturbing effects 
of the earth, sun, and solar radiation are modeled explicitly 
in this orbit determination process. Kepler element variations 
estimated by the process are ascribed to non-central 
lunar gravitational features. Doppler data are reduced in this 
manner for as many orbits as are available. In the second step, 
lunar gravity coefficients are determined using another weighted 
least-squares processor which fits the long-period Lagrange 
perturbation equations to the estimated Keplerian rates. 

Pseudo Doppler data have been generated simulating 
two different lunar trajectories. The perturbations included 
were triaxial lunar gravity harmonics, the earth, sun, and 
solar radiation pressure. Orbit determinations were performed 
using the empirical processor and the long-period orbital 
element variations obtained. The Kepler element rates from 
these convergences were used to recover the triaxial lunar 
gravity coefficients. The overall results of this controlled 
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experiment show that lunar gravity coefficients can be accurately 
determined and that the method as a whole is dynamically consist- 
ent. The pseudo data analysis shows the necessity of having 
Doppler data from different orbital inclinations for good 
selenodesy results. 

The method has been applied to Doppler data from the 
Lunar Orbiter I, II, III, and V missions. One hundred ninety- 
nine sets of Kepler element rates were obtained for lunar gravity 
field determination. A gravity field of degree and order four 
is derived from these Kepler element rates. Equipotential 
surfaces from this gravity field show the lunar mass distribution 
to be that of a triaxial ellipsoid with three large areas of mass 
concentration. The largest and by far the most dominant of 
these areas is centered very near the Mare Serenitatis region 
and covers a large portion of the near side of the moon. The 
other two regions of mass concentration are located on the far 
side of the moon but do not correspond to any specific mare 
region . 

This gravity field has been investigated using data 
from several of the Apollo missions. Orbit determination solu- 
tions (using a standard least-squares processor) from these data 
show that this fourth degree field results in improved orbit 
predictions as compared to those using other lunar gravity fields. 
All solutions indicate the lunar field models are still incomplete. 
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There are three major areas of this investigation 
which are original. The first is the development and applica- 
tion of an empirical orbit determination method for lunar orbits. 

The second is the derivation of a selenodesy method based on 
empirically determined Kepler element rates. The third is the 
generation of a fourth degree and order lunar gravity field 
(presented in this study) from Lunar Orbiter data using this 
method. This is the only indirect type selenodesy method that 
truly estimates long-period orbital element variations . 

This study demonstrates the necessity of obtaining 
more tracking data from lunar satellites at different orbital 
inclinations (specifically in the I = 26° to 85° gap). Only 
when this has been accomplished will correlations diminish 
between the various spherical harmonic coefficients and then 
an accurate gravity field will be determined. 

It is recommended that a version of this method 
which is determinate for both near-circular and near-equatorial orbits 
orbits be applied to Doppler tracking data from the Apollo 15 
Sub-satellite and any future lunar orbiters. 



APPENDIX A 




TRANSFORMATION FROM ORBITAL ELEMENTS TO STATE 


The rectangular state components (r,v) are found 
from the orbital elements, ic, using the following set of 
non-linear relationships : 


r = [R t ] g 


(A-l) 


where the entries to the R T matrix are : 


r ll 

= 

COS 

a 

cos 

(0 

- 

sin 

a 

cos 

I 

sin 

(0 

r l2 

= 

- cos 

ft sin 

to 

- sin 

a 

cos 

i 

COS to 

r 13 

- 

sin 

a 

sin 

I 










r 21 
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sin 
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cos 

0) 

+ 

cos 
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cos 
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sin 
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r 2 2 
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- sin 

ft sin 

to 

+ cos 
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cos 

I 

COS (0 

r 23 

= 

- cos 
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r 31 
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sin 

00 










r 32 

= 

sin 

i 

cos 

to 










r 33 

= 

cos 

i 













\ (A-2) 


and : 


q = 


a (cos E - e) 

I 2 . _ 

a/l-e sm E 


(A-3) 


0 
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The eccentric anomaly, E, used in these equations 
is found, given M and e, by solving Kepler's equation. 

E - e sin E = M (A-4) 


The velocity is found as follows: 

V = [R t 1 q 


(A-5) 


where : 


q = 


- sin E 


/l-e 2 


cos E 


V y « A 


(l-e cos E) 


(A-6) 


Hence the state x can be found from k at any time using 
A-l and A-5. 

TRANSFORMATION FROM STATE TO ORBITAL ELEMENTS 

Given (r,v) , the orbital elements are obtained as 

follows : 


h = r x v 


(A- 7 ) 


fi 



/ 




(A-8 ) 


I 


(A-9) 
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where h 


1 ' 


are the components of h. 


P = [S] r 


(A-10) 


where the entries to the S matrix are: 


then : 


S 11 

- 

COS 

n 





S 12 

= 

sin 

Q 





S 13 

=: 

0 






S 21 


- sin 

a 

cos 

I 

S 22 

= 

cos 

a 

cos 

I 


S 2 3 

— 

sin 

i 





S 31 

= 

sin 

n 

sin 

I 


S 32 

= 

- sin 

T 

cos 

n 

S„ 

33 

= 

cos 

i 






a) + f = tan 1 |p 2 / P^ 

r 2 h 2 / 2 i i/ ' 2 

r = v - h / r 
r ^ 

a “ 2 
(2y.- rv ) 

e = /l - h 2 / y^a 




; (a— ii ) 




(A-12) 


(A-13) 

(A-14) 


(A-15 ) 


cos 



- r) 
ae 


(A-16 ) 
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sin E 


rr 

e /y^a 


f 


tan 


■[ 


/l - & sin E 
cos E - e 


M = E - e sin E 


(A-17) 


(A-18) 


(A-19 ) 



APPENDIX B 


DOPPLER MODEL EQUATIONS 

The Doppler frequency shift data is modeled using 
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the range difference method. The estimated Doppler observable, 

p , can be computed over some counting interval (typically 
sixty seconds) r, using the following set of range equations 
(see Figure B-l) . 


P 4 = r* 


, ^r 


- r 


sr 


(t r ) 


(B-l) 



p 2 


= r* 



- r 

sr 


(t r - t) 


(B-2 ) 


(B-3 ) 


P 1 


r* 







1 

i 

i 

j 


(B-4) 


and 


r* = r + r 


em 


(B-5> 


where r is the earth-moon distance. p,, p„, p_, p. are 
em 12 K 3 4 

the topocentric ranges of the satellite, r sr and r . 


are the 
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receiving and transmitting station positions in earth centered 
inertial coordinates, t r is the Doppler time (at the end of 
the counting interval) and is the time the signal is acquired 
by ttfe receiving station, c is the speed of light, and t is 
the length of the counting interval. 

The equations for these four ranges are solved 
iteratively. The estimated observable is found as follows: 


1 


f , \ ( \ *1 

p = 27 

L 

p 3 + p 4) ‘ P 1 + p 2) 



\ / \ / J 


(B-6 ) 


where p is expressed in units of length per time. This value 

. . x p 4 

of p is considered valid at the time t = t - - — . 

r 2 c 



APPENDIX C 


PARTIAL DERIVATIVES 

I. PARTIAL DERIVATIVES OF DOPPLER WITH RESPECT TO STATE 

The partial derivatives of the Doppler data with 

40 

respect to satellite position and velocity are: 


3p_ _ 

3r 



(C-l) 


3p_ 

3 v 



(C-2 ) 


In these equations: 


P/i “ r ftr 2 c~ 


r sr ,tr " 


(C-3) 


1 1 * I 
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p 4 = r * ( - 2 


t p 4 
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- r_„ I t, - i 


sr 
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(C-4) 


P -3 = r* t 


_ t _ 4 

r 2 ~ c 




r st 1 ^ 2 


(C-5 ) 


p, = r* |tr - 2 - c 


- r 


st t tr 2 


x If 3 + p 4 j 


(C-6 ) 


dp. p. • p. 


and p ± = | Pi l, 


i a 
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II. 


(r,v) 


where 


PARTIAL DERIVATIVES OF STATE WITH RESPECT TO KEPLER 
ELEMENTS 


The partial derivatives of the satellite state 
with respect to the Kepler elements (a,e, 1,9 , u>,M) are: 


3r _ r 
3 a a 


(C-7) 


H- If 
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3e 


3r 
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r (1-e ) 

q iq 2 

r (1-e 2 ) 
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- z cos 
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X 

0 
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where : 
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= 
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t 23 
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= 
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9v 

81 


z sin B 
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y cos fi - x sin fi 


9 v 

9 Q 


- y 

X 

0 


(C-19 ) 
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where 



III. PARTIAL DERIVATIVES OF THE ORBITAL ELEMENTS WITH RESPECT 
TO THE KEPLERIAN PARAMETERS 

A typical orbital element, eccentricity for example, 
is represented in the following way: 


e(t) = e + e.t + e„t‘ 
o 1 2 


+ 6e 


Se o + 4e SR 


(C-23) 


Although the variations arising from the earth, sun. 


and solar radiation are functions of the Kepler state variables. 



137 


analysis has shown this dependence to be extremely weak. As 
a result the partial derivatives of the Kepler elements with 
respect to the solution parameters do not include the third 
body and solar radiation effects. It should be noted that the 
omission of these small parts (on the order of one ten thousandth 
the smallest existing term) does not affect convergence. Hence 
the necessary partial derivatives are: 
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APPENDIX D 


MEAN VALUE SEMI-MAJOR AXIS EQUATION FOR 
EARTH PERTURBATIONS 


The mean value calculation for the semi-major axis 

of a lunar satellite under earth perturbations is carried out 

41 

using a procedure following Danby ' s for the sun perturbations 
on the moon's orbit about the earth. Essentially this method 
finds an equivalent mass for the three body system and using 
this mass value and the Kepler motion equation, deduces an 
adjusted semi-major axis value (mean value) . The accuracy 
to which this derivation is valid is of order e. Hence all 
terms that are eccentricity dependent are assumed zero. This 
assumption introduces only small errors for the purposes of 
this investigation since all orbits used have an eccentricity 
of e = . 3 or less. 

The radial perturbing force experienced by the satellite, 
averaged over the satellite and earth periods is given as 
follows : 


2 

n a 


F = 
r r 


[1 - 3/2 sin I] 


2 


(D-l) 
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The total averaged radial force exerted on the 
lunar satellite is as follows: 


H - 

— — - F 
2 r 
a 


(D-2 ) 


If an equivalent mass, y, is defined this radial force can then 
be obtained from an equivalent central force law. Hence 


2 “ 2 r 

a a 


(D-3) 


or 


y = y i 


n 

1 + (3/2 sin I - 1) 

2n^ 


( D— 4 ) 


The mean value of the mean motion for third body perturbations 

42 

is given by Anderson and is as follows: 
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n = n 


L 


7 n © 2 

1 - T ~2 d - 3/2 sin I) 

n 


(D-5 ) 


where 


n ° 


and a is the average value (constant in this 


case) of the semi -major axis. Using equations D-4 and D-5 an 
equivalent Kepler motion law can be written for the perturbed 


motion. 
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-3 -2 'v 
a n = y 


(D- 


where a is the desired mean value of the semi-major axis, 
this equation using the expressions for n and y, the mean 
value expression for the semi-major axis is as follows: 
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a = a 


1 + 
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(1 - 3/2 sin I) 


(D- 


6 ) 

Solving 

7 ) 



APPENDIX E 


WEIGHTED LEAST-SQUARES DATA REDUCTION 


The covariance matrix for the Doppler data weighted 
least-squares process, assuming zero mean errors in the 
estimates, is found by forming the expected value of {^AK AK T } 
using equation (3.31); 


Exp { AK AK T }= (H T WH ] -1 H T W Exp {Ap Ap T }wH[H T WH] 1 (E-l) 


If it is now assumed that the residuals are serially uncorrelated 
and are normally distributed random variables such that 


Exp ^ Ap Ap T }= 


(E-2) 


2 ^ 
(where a is the variance of the Doppler measurements and I xs 

the identity matrix) and the weight used in the least-squares 

processing is the inverse variance of the Doppler, 


W = [a 2 I] 1 


(E-3) 


then 


Exp { AK AK T > = 


T -1 
[H WH] 


(E-4 ) 
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T -1 

It is only under these conditions that the [H WH] matrix is 
the true covariance matrix of the process. 

If conditions (E-2) and (E-3) are not satisfied, the 

T -l 

[H WH] matrix is not the covariance matrix. For the case of 

the empirical orbit determination processor neither of these 

conditions can be met since the residuals are systematic and 

are also serially correlated. Hence no interpretation of 

T -1 

variance is made of the diagonal terms in the [H WH] matrix. 



BIBLIOGRAPHY 


1. Lorell, J. , "Lunar Orbiter Gravity Analysis," The Moon, 

Vol. 1, 190-231, 1970. ‘ " ' 

2. Michael, W. H. , "Physical Properties of the Moon as 

Determined From Lunar Orbiter Data," Presented at the 
Fourteenth General Assembly of the International Union 
of Geodesy and Geophysics Meeting, Lucerne, Switzerland, 
Sept. 1967. 

3. Tolson, R. H. and Gapcynski, J. P. , "An Analysis of the 

Lunar Gravitational Field as Obtained from Lunar 
Orbiter Tracking Data," Presented at the IQSY/COSPAR 
Assemblies, London, England, July, 1967. 

4. Gapcynski, J. P. , Blackshear, T. W. and Compton, H. R. , 

"The Lunar Gravitational Field as Determined from the 
Tracking Data of the Lunar Orbiter Series of Spacecraft," 
Presented at the AIAA-AAS Astrodynamics Specialists 
Conference, Grand Teton National Park, Wyoming, Sept. 1968. 

5. Blackshear, T. W. , Compton, H. R. , and Schiess , J. R. , 

"Preliminary Results on the Lunar Gravitational Field 
from Analysis of Long-period and Secular Effects of 
Lunar Orbiter I," Presented at NASA Seminar on Guidance 
Theory and Trajectory Analysis, Electronics Research 
Center, Cambridge, Mass., May 1967. 

6. Lorell, J. and Sjogren, W. L. , "Lunar Gravity: Preliminary 

Estimates from Lunar Orbiter," Science, Vol . 159 , 3815 , 1968. 

7. Risdal, R. E., "Development of a Simple Lunar Model for 

Apollo," Contract Report D2-100819-1, The Boeing Company, 
Seattle, Washington, 1968. 

8. McCuskey, S. W. , Introduction to Celestial Mechanics , Addison- 

Wesley Publishing Company, Inc . , Reading , Mass . , 1963. 

9. Caputo, M. , The Gravity Field of the Earth from Classical 

and Modern Methods , Academic Press, New York, 1967. 

10. "Transactions of the International Astronautical Union," 

Vol. XIB , Proceedings of the 11th General Assembly , 

Berkeley, 1961, D« H. Sadler ed.. Academic Press, New York, 
1962 . 


144 


11. Kaula, W. M. , Theory of Satellite Geodesy , Blaisdell 

Publishing Co., Waltham, Mass., 1966. 

12. Pollard, H. , Mathematical Introduction to Celestial Mechanics , 

^ Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1966. 

13. Bullock, M. V. and Ferrari, A. J. , "Orbit Determination 

For Lunar Parking Orbits Using Time-Varying Orbital 
Elements," NASA Contractor Report 110008, May 1970. 

14. Kozai, Y. , " The Motion of a Close Earth Satellite," Astro. J., 

Vol. 64, 378-397, 1959. 

15. Kozai, Y. , "Second Order Solution of Artificial Satellite 

Theory Without Drag." Astro. J . , Vol. 67, 446- 61, 1962. 

16. Brouwer, D. and Clemence , G. M. , Methods of Celestial 

Mechanics , Academic Press, New York, 1961. 

17. Lorell, J. and Liu, A., "Method of Averages Expansion for 

Artificial Satellite Application," Jet Propulsion Lab. 
Report -32-1513, April 1971. 

18. Ferrari, A. J. and Heffron, W. G., "Effects of Physical 

Librations of the Moon on the Orbital Elements of a 
Lunar Satellite," Presented at AIAA-AAS Astrodynamics 
Conference, Ft. Lauderdale, Fla., August 1971. 

19. Kozai, Y. , "Effects of Solar Radiation Pressure on the 

Motion of an Artificial Satellite," Smithsonian 
Institution, Special Report No. 56, January 1961. 

20. Gaposchkin, E. M. , "Differential Orbit Improvement," 

Smithsonian Institution, Special Report No. 161, 1964. 

21. de Vezin, H. G. , "Doppler Observable Modeling for Apollo 

Real-Time Orbit Determination Program," Presented at 
Astrodynamics Conference, Manned Spacecraft Center, 

Houston, Texas, Dec. 1967. 

22. Sage, A. P. , Optimum Systems Control , Prentice-Hall, Inc., 

Englewood Cliffs, New Jersey, 1968. 

23. Izsak, I. G. , "Differential Orbit Improvement with the 

Use of Rotated Residuals," Smithsonian Institution, 

Special Report No. 73, 1961. 

24. Wollenhaupt, W. R. , "Apollo Orbit Determination and 

Navigation," Presented at AIAA 8th Aerospace Sciences 
Meeting, New York, N. Y. ,- January 1970. 



145 


25. Kozai, Y. , op. cit., 1959. 

26. Kozai, Y. , Smithsonian Astrophysical Observatory, Private 

Communication, February 1971. 

27. Anderson, J. D. , "Long Term Perturbations of a Moon Satellite 

by the Earth and Sun," Jet Propulsion Laboratory 
Technical Memorandum 312-162, February 1962. 

28. Durbin, J. , "The Fitting of Time - Series Models," Revue 

Inst. Int. De Stat . , 233-243, 1960. 

29. Battin, R. H. , As tronautical Guidance , McGraw-Hill Book Co., 

New York, 1964. 

30. Melbourne, W. G. et al, "Constants and Related Information 

for Astrodynamic Calculations, 1968," Jet Propulsion Lab. 
Technical Report 32-1306, July 1968. 

31. "Natural Environment and Physical Standards for the 

Apollo Program and .the Apollo Applications Program," 

NASA M-DE-8020 . 00C , SE -15-001-lB, July 1969. 

32. Peabody, P. R. , Scott, J. F. , and Orozco, E. G. /'Users 

Description of JPL Ephemeris Tapes," Jet Propulsion 
Laboratory Report No. 32-580, March, 1964. 

33. Muller, P. M. , and Sjogren, W. L. , "Lunar Mass Concentrations," 

Science , Vol. 161, 3842, 1968. 

34. Kaula, W. M. , op. cit. , 1966. 


35. Koziel, K. , "The Constants of the Moon's Physical Libration 

Derived on the Basis of Four Series of Heliometric 
Observations from the Years 1877 to 1915," Icarus, 

Vol. 7, 1-28, 1967. 

36 . Michael, W. M. , et al , "Results on the Mass and Gravitational 

Field of the Moon as Determined from Dynamics of Lunar 
Satellites," Dynamics of Satellites 1969 , Bruno Morando 
ed. , Springer-Verlag , Berlin, 1970. 

37. Koziel, K. , op. cit. , 1967. 


38. Liu, A. S. and Laing , P. A., "Lunar Gravity Field as Deter- 
mined by Orbiters , " Presented at 14th Plenary Meeting 
of COSPAR, Seattle, Washington, June 1971. 


146 


39. 

de Vezin, H. 

G. ' 2E* cit. , 1967. 

40. 

de Vezin, H. 

G. , o£. cit. , 1967. 

41. 

Danby, J. M. 
Macmillan 

A., Fundamentals of Celestial Mechanics, 
Company, New York, 1962. 

42. 

Anderson, J. 

D. , o£. cit. , 1962. 



